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A  GENERALIZATION  OF  THE  LR  ALGORITHM 
TO  SOLVE  AX  =  \  EX 

BY 

Linda  Kaufman 

Abstract 

In  this  paper,  we  will  present  and  analyze  an  algorithm  for 
finding  x  and  \  such  that 

<-w» 

Ax  =  xBx  (1) 

^  /X/ 

where  A  and  B  are  n  x  n  matrices.  The  algorithm  does  not  require 
matrix  inversion,  and  may  be  used  when  either  or  both  matrices  are 
singular.  Our  method  is  a  generalization  of  Rutishauser's  LR  method  [17] 
for  the  standard  eigenvalue  problem  Ax  =  \x  and  closely  resembles  the 

^  l*W 

QZ  algorithm  given  by  Moler  and  Stewart  [10]  for  the  generalized  problem 
given  above.  Unlike  the  QZ  algorithm,  which  uses  orthogonal  transfor¬ 
mations,  our  method,  the  LZ  algorithm,  uses  elementary  transformations. 
When  either  A  or  B  is  complex,  our  method  should  be  more  efficient. 
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or  in  part  is  permitted  for  any  purpose  of  the  United  States  Government. 


The  LZ  algorithm  is  based  on  three  observations: 

1)  If  L  and  M  are  nonsingular  matrices,  the  eigenvalue  problems 
LAMy  =  XLBMy  and  Ax  =  XBx  have  the  same  eigenvalues  and  their  eigen¬ 
vectors  are  related  by  x  =  My. 

2)  If  A  is  a  triangular  matrix  with  diagonal  elements  cr* , 
and  B  is  a  triangular  matrix  with  diagonal  elements  p^,  then  for  each 

3,  i  =  l,2,---,n,  is  an  eigenvalue  of  the  generalized  problem  if  p^O. 

If  for  some  i,  is  zero,  then  the  polynomial,  determinant  (A-Xb)  > 
is  of  degree  less  than  n.  If  is  not  zero  and  the  corresponding 
p^  is  zero,  we  say  that  infinity  is  an  eigenvalue.  If  for  i} 

both  <y^  and  p^  are  zero,  then  det(A-XB)  vanishes  for  all  values 
of  X,  and  every  scalar  is  an  eigenvalue  of  Ax  =  XBx  . 

3)  There  exist  matrices  I.  and  M  such  that  LAM  and  LEM  are 
upper  triangular  and  L  and  M  are  the  products  of  lower  triangular  and 
permutation  matrices. 

The  first  two  observations  should  be  obvious;  the  third  requires 
explanation.  In  [18]  Stewart  shows  that  there  exist  two  unitary  matrices 
U  and  V  such  that 

A'  =  l^AV  and  B/  =  l^BV 

are  upper  triangular.  The  symbol  indicates  the  conjugate  transpose 
of  the  matrix  U.  We  can  certainly  write 

as  RL  and  V  as  MS 

vhere  S  and  R  are  both  upper  triangular  matrices  and  L  and  M  are 
products  of  lower  triangular  and  permutation  matrices.  The  matrices 
R  ^A 'S  ^  =  LAM  and  R  ^B'S  ^  =  LBM  ere  both  upper  triangular  and  hence 


verify  our  observation. 


t 


The  LZ  algorithm  has  three  parts.  In  the  first  part,  the  matrix 

A  is  reduced  to  upper  Hessenberg  form,  i.e.  a. .  =  0  for  i  >  j+L, 

■*■  3 

and  B  is  simultaneously  reduced  to  triangular  form.  The  first  stage 
of  the  LZ  algorithm  is  very  similar  to  the  first  stage  of  the  Moler- 
Stewart  algorithm,  and  they  may  be  freely  substituted  for  each  other. 

The  advantage  of  using  our  method  is  that  it  is  about  5/2  faster;  the 

advantage  of  theirs  i3  numerical  stability.  The  second  stage  of  the  LZ  algorithm 

is  a  generalization  of  the  LR  algorithm  and  iteratively  reduces  A  to 

triangular  form  while  preserving  the  triangularity  of  B.  The  last  part 

of  LZ  obtains  the  eigenvectors  of  the  triangular  matrices  and  transforms 

them  back  into  the  original  coordinate  system.  Throughout  the  algorithm 

stabilized  elementary  transformations  (see  Wilkinson  [ 19] ,  p.  16U)  are 

used  to  insure  numerical  stability.  These  transformations  are  the 

products  of  lower  triangular  matrices  and  permutation  matrices,  and 

are  easy  to  compute  and  easy  to  use.  The  permutation  matrices  are 

designed  to  help  minimize  the  loss  of  accuracy  in  numerical  operations. 

A  further  explanation  of  the  stabilized  elementary  transformations  used 
in  the  heart  of  the  LZ  algorithm  is  contained  in  the  notation  section 
at  the  end  of  this  introduction. 

I.  BACKGROUND 

As  Lancaster  (  8]  and  Gantmacher  [6  ]  point  out,  the  generalized 
eigenvalue  problem  often  occurs  in  the  physical  sciences.  Many  mechan¬ 
ical  and  electrical  systems  are  governed  by  a  differential  equation  of 
the  fora 

Cx  +  Dx  +  Ex  =  0 


a 

i 

i 


Ji 


9 
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where  C,  D,  and  E  are  n  x  n  matrices  and  the  solution  is  expected  to 
haw  the  form  x(t)  =  e^Sc(O) .  Solving  the  ordinary  differential  equation 
entails  finding  the  eigensystem  of 

(X2C  +  XD  +  E)x  =  0 


If  no  damping  occurs  and  the  D  term  is  missing,  the  problem  is  like 

p 

the  one  given  in  (l)  in  X’-  If  the  system  is  damped  and  the  D  term 
is  present,  the  problem  can  be  reconstructed  to  have  the  form  of  (1) 
where  now  A  is  the  matrix 


and  B  is  the  matrix 


In  many  problems,  especially  those  which  describe  physical  systems, 

A  and  B  have  some  special  structure  and  most  of  the  algorithms  in  the 
literature  are  designed  for  matrices  having  specific  properties.  In  [  9* 
Martin  and  Wilkinson  have  given  a  method  for  A,B  symmetric  and  B  positive 
definite.  Crawford  [  2)  nee  presented  a  modification  of  that  algorithm 
when  B  is  a  band  matrix.  In  [ 7 ]j  Golub,  Underwood  and  Wilkinson  describe 
a  version  of  the  Lanczos  algorithm  for  A,B  symmetric  and  B  positive  definite. 
Fix  and  Heiberger  [ 3 ]  have  a  method  designed  for  illconditioned  B  which 
requires  the  determination  of  the  rank  of  certain  submatrices  in  A  and  B. 

If  symmetry  and  positive  definiteness  are  not  present  and  B  is  well  con¬ 
ditioned,  the  eigensystem  of  Ax  =  XBx  can  be  found  by  forming  B  and 
determining  the  eigensystem  of  3  ^Ax  =  Xx,  for  wnich  there  exist 
several  good  methods.  For  a  nearly  singular  B,  Peters  and  Wilkinson  [15] 


U 


describe  an  algorithm  which  approximates  the  null  space  of  B.  This 


approach  involves  determining  the  rank  of  submatrices  which  is  often  difficult 
to  do  exactly  on  a  finite  precision  computer. 

Recently  Moler  and  Stewart  [  10]  have  presented  an  algorithm  for 
solving  the  generalized  eigenvalue  problem  which  may  be  used  regardless 
of  the  cornition  and  structure  of  the  two  matrices.  Our  algorithm 
resembles  the  Q£  algorithm  in  that  we  generalize  Rutishauser 's  LR 
algorithm  [  17]  in  the  same  way  that  Moler  and  Stewart  generalize  Francis's 
QR  method  [ 5  ]  for  the  standard  eigenvalue  problem  Cx  =  Xx.  Before 
we  describe  our  algorithm  in  detail  and  discuss  its  relationship  to  the 
Q2  method  perhaps  it  is  best  to  review  the  Qft  and  the  LR  methods.  In 
practice,  the  LR  method  for  the  problem  Cx  =  Xx  is  essentially: 

1)  Reduce  C  to  upper  Hessenberg  form  using  similarity 
transformations . 

2)  Determine  a  shift  X- 

3)  Find  L,  a  product  of  stabilized  elementary  transfor¬ 
mations,  and  R,  an  upper  triangular  matrix,  so  that 
L(C  -  XI)  =  R. 

M  Set  C  to  LCL-1.  The  matrix  C  will  oe  upper  Hessenberg. 

5)  If  the  subdiagor.al  elements  of  C  are  r-.ot  negligible. 


return  to  2 . 

6)  The  eigenvalues  of  the  original  matrix  are  the  diagonal 
elements  of  C. 
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according  to  the  ratios  of  the  shifted  eigenvalues.  In  practice,  the  shift  i 
usually  an  eigenvalue  of  the  lowest  2x2  subblock  on  tne  diagonal  of 
C  which  has  not  been  t riangu lari zed .  This  policy  often  gives  a  good 


approximation  to  an  eigenvalue  of  the  whole  matrix. 

The  basic  QR  method  is  approximately  the  same  as  the  LR  method 
with  an  or^ogonal  matrix  Q  replacing  the  matrix  L  in  steps  3  and  4. 
In  practice  a  double  shift  implicit  version  of  the  QR  method  is  used  in 
which  steps  3  and  4  read 

3)  Find  Q,  an  orthogonal  matrix,  and  R,  an  upper  triangular 

matrix,  such  that 

R  =  QT  =  Q(C  -  XI)  (C  -  XI)  where  X  and  X 

are  complex  conjugate  shifts  or  a  pair  of  real  shifts. 

h)  Set  C  to  QCQT. 

Only  the  first  column  of  T  is  ever  explicitly  formed. 

The  main  advantage  of  the  double  shift  algorithm  is  the  preser¬ 
vation  of  real  arithmetic  for  real  matrices.  The  QZ  algorithm  also  has 
this  property.  With  the  double  step  QJR  and  QZ  methods  the  final  matrix 
is  not  necessarily  triangular,  but  may  have  2x2  submatrices  on  the 
diagonal  which  must  be  resolved.  The  LR  and  the  LZ  algorithms  do  not 
limit  themselves  to  real  arithmetic  but  avoid  the  2x2  problem.  Double 
shift  LR  and  LZ  methods  are  not  found  in  practice  because  of  the  lack  of 
a  theoretical  basis.  Francis  [5  ]  has  proved  that  one  iteration  of  the 
implicit  double  shift  QR  method  is  equivalent  to  two  iterations  of  the 
basic  QR  method.  His  theorem  is  based  on  the  uniqueness  of  orthogonal 
transformations  which  reduce  a  given  matrix  to  triangular  form  with 
positive  diagonal  elements.  This  uniqueness  property  is  missing  for 
stabilized  elementary  transformations. 
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The  LZ  Algorithm  as  a  Generalization,  of  the  LR  Algorithm 


The  LZ  algorithm  is  motivated  by  the  LR  method  described  above  where 
the  matrix  C  is  AB  However,  we  do  not  assume  that  B  ^  exists. 


Briefly,  our  algorithm  is: 


1)  Reduce  A  to  upper  Hessenberg  form  and  B  to  triangular 

form.  : 

2 )  Find  a  shift  X • 

3)  Find  matrices  L  and  •  M,  stabilized  elementary  trans¬ 
formations,  such  that  L(A  -  XB)  is  upper  triangular  and 
LBM  is  upper  triangular. 

U)  Set  A  to  LAM  and  B  to  LBM.  The  new  A  will  be  upper 
Hessenberg. 

5)  If  the  subdiagonal  elements  of  A  are  not  negligible,  return  to  2 

6)  The  i  1  eigenvalue  is  a../b..  if  b. .  is  nonzero. 

ir  11  11 

Again  the  shift  is  used  to  hasten  che  convergence  of  the  algo¬ 
rithm.  In  practice  it  is  usually  a  solution  of  the  lowest  2x2  sub¬ 
problem  on  the  diagonal  of  A  -  XB  which  has  not  been  triangularized. 

If  the  matrix  C  in  the  LR  method  is  AB  1,  then  the  matrix  L 
in  the  third  step  of  the  LR  method  is  precisely  the  matrix  L  in  the 
third  step  of  the  LZ  method  if  both  algorithms  employ  the  same  pivoting 
strategy.  This  fact  is  verified  by  denoting  the  left  hand  transformation 
in  the  third  step  of  the  LZ  algorithm  by  L  and  noticing  that 


L(C  -  XI)  =  L(AB  -  XI) 
=  L(A  -  XB)B_1 


an  upper  triangular  matrix.  Thus  L  is  also  a  transformation  which 


J 


triangularises  C  -  \I;  and  if  the  pivoting  strategy  is  the  same,  L 
is  the  transformation  L  in  step  3  of  the  method  LR. 

Moreover,  it  can  be  shown  that  che  two  algorithms  produce 
corresponding  iterates.  If  C'  denotes  the.  next . iterate  in  the  LP. 
algorithm  and  A'  and  B'  are  the  iterates  in  the  LZ  algorithm, 
then  in  the  LR  algorithm  ■ . 

C'  =  lcl"1 

and  in  the  LZ  algorithm 

A'b'"1  =  LAM4‘1B“1L"1 
=  LAB-iL-1 
=  LCL-^ 

•  I 

In  Chapter  One  we  will  present  «wo  algorithms.  The  first  is  a 
straight  generalization  of  the  LR  method.  The  second  is  an  implicit 
scheme  in  which  only  the  first  column  Of  A  -  lB  is  actually  formed. 
The  second  scheme  requires  fewer  operations  and  is  more  stable. 

III.  NOTATION 


To  simplify  the  explanation  in  the  remainder  of  this  paper,  we 

introduce  the  following  symbols ; 

For  a  complex  scalar  a,  lid’ll  will  denote  jlra(d)|  +  |Re(d)|  . 

||d||  corresponds  to  the  1  norm  of  ot,  if  d  is  considered  as  a  vector 
in  the  complex  plane. 


In  general,  the  (i,j)th  element  of  the  matrix  A  will  be  denoted 


by  a^.  If  a  matrix  A  is  the  kth  element  in  a  sequence  of  matrices, 
it  will  be  designated  by  A and  its  (i,j)  element  will  be  designated 
by  . 
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£^  will  denote  the  subset  of  the  set  of  stabilized  elementary 
matrices  described  in  Wilkinson  [  19]  having  the  form 


t 


t 


» 


t 


t 


where  ||t^|!  <  1.  Blank  spaces  indicate  zeroes. 

When  a  matrix  is  multiplied  on  the  left  by  a  matrix  of  £.  of  the  first  form 

J. 

s  t 

only  its  i+1  row  is  changed,  but  when  a  matrix  is  multiplied  on  the 

left,  by  a  matrix  of  £  of  the  second  form,  its  ith  and  i+lSt  rows  are  first 

interchanged  and  then  a  raultip  .e  of  the  new  itn  row  is  added  to  the 
•  1 st 

i.+l  row. 


We  will  often  use  a  member  of  £^  to  annihilate  an  element  in 

S 

the  i+1  row  of  a  matrix.  For  example,  we  may  want  to  zero  a.  ,  .  . 

If  either  the  current  a.  .  or  a. .  is  nonzero,  then  there  exists 

i+J-  >0  1  j 

a  unique  member  of  £_.  which  will  annihilate  a4il  .  .  Specifically, 

if  j|  ai+^  j|  is  less  than  or  equal  to  y  a^  j|  ,  we  use  a  matrix  of  the  first 

form  where  !).  is  -a..,  ./a.  .  ;  otherwise,  we  use  a  matrix  of  the 

1  i+J-,0  ij 

second  form  with  T).  given  by  -a.  ./a...  .  .  If  both  the  current  a.  ,  . 

x  ij  l-rx,j  x^Xmj 

and  a  are  zero,  then  any  member  of  £.  will  leave  a  zero  in  a.  .  .  ■ 


% 

form 


v.here  |j^u^jj  <  1,  or  the  form 

.th 

1  row— 

where  U^J  <  i. 

When  a  matrix  is  multiplied  on  the  right  by  a  matrix  of  the  first 

form,  only  its  i  column  is  changed,  but  when  a  matrix  is  multiplied 

t/h  s"t 

on  the  right  by  a  matrix  of  the  second  form,  its  i  and  i+1 

s  t 

columns  arc  interchanged  and  a  multiple  of  the  new  i+1  column  is 
added  to  the  j  '*  column. 

The  set  T  will  denote  the  set  of  matrices  in  upper  triangular 

form.  If  A  is  in  T  ,  then  a.  .  =  0  for  i  >  j.  The  set  V  *111 

1  jJ 

denote  the  set  of  matrices  in  uppei-  Hessenberg  for~  If  A  is  in  y  , 

then  a.  .  =  0  for  i  >  j+1. 
id 

Each  iteration  of  the  LZ  algorithm  invj_vres  multiplying  matrices 
by  product  of  transformations.  In  our  discussion  of  the  LZ  algorithm 
the  symbol  A '  will  usually  denote  the  matrix  A  after  all  the  trans¬ 
formations  for  one  iteration  have  been  applied  to  it.  The  symbol  A* 
will  represent  the  matrix  A  after  some  but  not  all  of  the  transfor¬ 
mations  for  one  iteration  have  been  applied  to  it. 


CHAPTER  ONE 


In  this  chapter  we  shall  describe  the  LZ  algorithm  in  detail. 

As  mentioned  in  the  introduction,  the  algorithm  has  3  sections: 

1)  Reducing  B  to  triangular  form  and  transforming  A  to  upper 
Hessenberg  form. 

2)  Iteratively  reducing  A  to  triangular  form  while  preserving 
the  triangularity  of  B. 

3)  Finding  the  eigenvectors  of  the  triangular  system  and 
transforming  them  into  the  eigenvectors  of  the  original 
system. 

To  obtain  the  eigenvectors  of  the  original  system  all  the  right 
hand  transformations  must  be  accumulated,  for  if  L  and  M  are  nonsingular 
matrices  and 

LAMjr  =  XLEMy 

then 

Ax  =  XBx 

where  x  =  My.  Thus,  if  y  is  an  eigenvector  of  the  triangular  system, 
is  an  eigenvector  of  the  original  system. 

In  the  second  section  of  this  chapter,  where  we  present  the 
iterative  section  of  the  algorithm,  we  will  describe  two  algorithms. 

The  first  method  is  an  explicit  scheme  which,  if  B  ^  existed,  would  be 
quite  like  the  LR  algorithm  for  AB  \  The  second  algorithm  is  a  more 

stable  implicit  scheme.  In  the  third  section  we  prove  that  the  two 
algorithms  generate  and  use  the  same  transformations. 


fc 
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a 


M 
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I.  INITIAL  REDUCTION  TO  HESSEN BERG -TRIANGULAR  FORM 

In  this  section  an  algorithm  will  be  described  that  reduces  a 
matrix  A  to  upper  Hessenberg  form  and  reduces  a  matrix  B  to  trian¬ 
gular  form  by  applying  the  same  elementary  transformations  to  both 
matrices. 

The  first  step  is  the  standard  Gaussian  elimination  process  with 
partial  pivoting  as  described  in  Forsythe  and  Mol-er-^  .  We  find  a 
matrix  L,  the  product  of  elementary  and  permutation  matrices  such  that 
LB  is  upper  triangular,  and  then  replace  A  and  B  by  LA  and  LB,  respec¬ 
tively. 

In  the  next  stage  we  reduce  A  to  an  element  of  V  while  main¬ 
taining  the  triangularity  of  B.  We  begin  by  choosing  an  element  Lr_^ 
from  £n  1  so  that  replacing  A  by  L^A  puts  a  zero  in  the  (n,l) 
position  of  A.  Multiplying  B  on  the  left  by  L  ^  introduces  a  new 
nonzero  element  in  the  (n,n-l)  position  of  B.  If  pivoting  had  been 
necessary,  B  would  still  have  the  same  form.  Thus  A  and  B  now  look 
like 

A  B 


X  X  X  X  X 

X  X  X  X  X 

X  X  X  X  X 

0  X  X  X  X 

X  X  X  X  A 

0  C  XXX 

X  X  X  X  X 

0  0  0  X  X 

0  X  X  X  X 

0  0  0  X  X 

We  now  focus  on  B,  and  choose  a  matrix  M  ,  from  n  so  that 
7  n-l  n-x 

setting  B  to  BM  .  and  A  to  AM  ,  returns  B  to  triangular  form  and  maintains 
n-i  n-J. 

the  zero  we  introduced  into  A.  Thus  we  have 
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e 
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A 

B 

X  X  X  X  X 

pc  X  X  X  X 

X  X  X  X  X 

oxxx'x 

i 

X  X  X  X  X 

0  0  X  X  X 

X  X  X  X  X 

•  0  0  0  X  X 

0  X  X  X  X 

0  0  0  0  X 

The  process  is  continued  until  A  is  in  V.  As  edch  element  in  '  , 

•  I 

A  is  zeroed  using  a  row  transformation,  a  nonzero  element  is  introduced 

t 

on  the  subdiagonal  of  B  vhi ch  must  be  immediately  annihilated  using  a 

•  l 

column  transformation.  Elements  are  eliminated :from  A  in  the  order  given  below 

X  X  X  X  X  ’ 

1  : 

X  X-  X  X  X  . 

X5-  X  X  X  X 

. '  X2  X5  X  X!  X  !  , 

XI  X4  X6  X  X 

j  1 

Note  that  pivoting  maintains  stability  but  does  not  affect  the 
zero  structure  of  the  two  matrices  any  more  than  a  'nonpivoting  algorithm 
would.  ‘ 

There  are  other  ways  to  annihilate  elements  of  A  and  B  which 
might  be  more  efficient  or  more  stable  for  any  given  problem.  One  such 
method  involves  reducing  to  B  to  an  element  of  T  and  then  using  column 
transformations  to  reduce  A  to  an  element  of  V.  The  nonzero  elements, 
which  are  introduced  on  the  subdiagonal  of  B  by  the  column  transfor¬ 
mations,  arc  eliminated  using  rew  transformations.  Elements  of  A  would 
be  zeroed  in  the  order  given  below: 
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X  X  X  X  X 


X  X  X  X  X 

:  6 

X  X  X  X  X 

:  4  5  - 

X  X^  X  X  X 

x1  x2  yp  x  x  . 

This  algorithm  would  certainly  be  more  efficient  than  the  first 
method  described  if  B  were  the  identity  matrix  and  if  A  were  the 
matrix 


■  /  1  1  1  1  1 
/  11110 

■  I  1  l'  1  0  0 

l  1  1  .0  0  0 
\1  0  0  0  0 

Both  algorithms  just,  described  require  about  13n3 */6  multipli¬ 
cations  and  additions.  In  terms  of  the  first  method  the  opera¬ 

tion  count  can  be  broken  down  in  the  following  way: 


1)  Reducing  B  to  triangular  form 

Transformations  on  A 
Transformations  on  B 

l 

2)  Reducing  A  to  an  element  of  V  and  pre¬ 
serving  the  triangularity  of  B 

a)  Tti  eliminate  elements  in  the  .i 
column  of  A 

Transformations  on  A 


Additions  +  Multiplications 

n5/2 

n3/3 


(2n-j)  (n-i-j) 
(n+2)  (n-l-j) 
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Transformations  on  B 


t  3 

\  i 


b)  Total  2nd  step 

Transformations  on  A  5n  /6 

Trans formal  ions  on  B  r^/2 

If  the  eigenvectors  are  requested,  then,  as  we  explained  at  the  begin¬ 
ning  of  this  chapter,  the  product  of  the  M’s  must  be  accumulated.  This 

requires  n5/2  additions  and  n5/2  multiplications. 

In  comparison,  the  first  part  of  the  QZ  algorithm  requires 

5  5  2 

17 nJ / $  multiplications,  17u/j  additions,  and  n  square  roots.  If 

eigenvectors  are  also  desired  .  the  QZ  algorithm  expends  an  additional 

2  2 
3n/2  multiplications  and  $nJ/2  additions. 

It  is  interesting  to  compare  the  above  figures  with  the  number 

of  operations  needed  to  fora  AB  1  and  reduce  this  matrix  to  Hessenberg 

form.  If  iterative  refinement  is  not  done,  the  basic  process  requires 

about  13 n?/6  multiplications  and  additions  or  the  same  number 

required  for  the  first  section  of  LZ.  If  eigenvectors  are  desired, 

x? / 2  extra  multiplications  and  x? / 2  additions  are  needed.  The  figures 

in  this  paragraph  assume  nonunitary  transformations  are  being  used. 

The  following  table  summarizes  the  cost  of  using  the  initial 

part  of  the  three  algorithms. 

summary  or  uperation  Counts 


Without  Eigenvectors 


With  Eigenvectors 


+ 

X 

square 

roots 

+ 

X 

square 

roots 

u 

| 

LZ 

l*>x?/6 
"  / 

13  n^/  6 

0 

l6n^/ 6 

iSx?/ 6 

0 

{ 

i 

r 

QZ 

$kx? /6 

34n5/6 

2 

n 

k^x? /S 

43n^/6 

2 

n 

t 

1 

AB_1 

13^/6 

l?n3/6 

0 

l6n^ /6 

l6n‘V  6 

0 

.  3 

*  > 

& 

II .  FINDING  THE  EIGENVALUES 


In  this  section  we  give  an  algorithm  for  determining  the  eigen¬ 
values  of  the  problem  Ax  =  XBx  where  A  is  upper  Hessenberg  and  B  is 
upper  triangular.  As  in  Moler's  and  Stewart's  Q£  algorithm,  the  method 
entails  iteratively  reducing  A  to  upper  triangular  form  while  preserving 
the  triangularity  of  B. 

Each  iteration  of  the  iterative  procedure  is  essentially: 

1)  Find  a  shift  X  which  could  be  a  _/b  or  an  eigen- 

nnr  nn 

value  of  the  lowest  2  by  2  subproblem  of  A  -  XB  . 

2)  Find  a  matrix  L  such  that  L(A  -  XB)  is  upper  triangular. 

5)  Find  a  matrix  M  such  that  LBM  is  upper  triangular. 

4)  Set  A'  to  LAM  and  B'  to  LBM.  A'  will  be  in  V. 

Most  of  this  sectic  i  discusses  the  construction  of  L  and  M  and 
their  application  to  our  matrices  to  satisfy  the  requirements  given  above. 
If  the  matrix  T  denotes  A  -  XB,  then  it  is  obvious  that  LAM  =  LTM  +  XLffi 
and  that  LAM  is  in  V  if  and  only  if  LTM  is  in  V. 

Each  iteration  begins  vith  an  A  which  has  r.o  zero  subdiagonal 
elements.  If  after  the  iteration  A'  has  a  zero  on  its  subdiagonal,  we  can 
deflate  the  problem  and  wo;k.  on  a  lower  dimensional  subproblem.  Hence 
the  purpose  of  each  iteration  is  to  drive  the  elements  on  the  subdiagonal 
A '  closer  to  zero.  In  Chapter  2  we  will  specify  conditions  under  which 
the  process, ve  are  about  to  describe, accomplishes  this  goal. 
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In  the  'explicit'  version  ol  the  algorithm,  we  start  an  iteration 
by  forming  T  ~  A  -  XB.  Since  A  is  upper  Hessenberg  and  B  is  upper 
triangulai ,  our  matrices  look  like 


X  X  X  X  X 


X  X  X  X  X 


X  X  X  X  X 


0  X  X  X  X 


0  X  X  X  X 


0  0  X  X  X 


0  0  X  X  X 


0  0  0  X  X 


0  0  0  X  X 


0  0  0  0  x 


We  now  select  an  element  from  I1  so  that  replacing  T  by  L^T  zeroes  the 
element  in  the  (2,1)  position  of  T.  Then  replacing  B  by  L^B  introduces  a 
nonzero  element  in  the  (2,1)  position  of  B.  T  and  B  now  have  the  form 


X  X  X  X  X 


X  X  X  X  X 


0  X  X  X  X 


X  X  X  X  X 


0  X  X  X  X 


0  0  X  X  X 


0  0  X  X  X 


0  0  0  X  X 


0  0  0  X  X 


0  0  0  0  x 


Our  attention  is  now  turned  to  3,  and  is  chosen  from  ^  so 
that  replacing  B  by  BM^  yields  a  triangular  matrix .  The  application  of 
to  T  is  delayed.  We  now  return  to  T  and  annihilate  'Ug  by  an 
element  1^  in  Applying  the  same  transformation  to  B  produces 


X  X  X  X  X 


X  X  X  X  X 


0  X  X  X  X 


0  X  X  X  X 


0  0  X  X  X 


0  X  X  X  X 


0  0  X  X  X 


0  0  0  X  X 


The  matrix  is  now  applied  to  T  and  is  chosen  from 
so  that  is  in  T.  The  matrices  T  and  B  now  look  like 


X  X  X  X  X 


X  X  X  X  X 


X  X  X  X  X 


0  X  X  X  X 


0  0  X  X  X 


0  0  X  X  X 


0  0  X  X  X 


0  0  0  X  X 


0  0  0  X  X 


0  0  0  0  X 


It  is  important  to  notice  that  the  element  t,n  is  zero. 

51 

Future  transformations  will  not  touch  this  element,  and  hence  it  will 
remain  zero  throughout  the  iteration,  furthermore,  the  situation  had 
not  been  influenced  by  the  form  of  M^,  i.e.  whether  pivoting  had  been 
necessary  to  stably  zero  b21  . 

In  general,  row  transformations  are  applied  to  T  and  B  simul¬ 
taneously  but  column  transformations  are  applied  first  to  B.  If  we 
write  M  as  .  .  .  Mn_1,  then  as  we  apply  to  B  we  apply  Mi_1  to 

T.  Each  row  transformation  will  zero  an  element  of  T  and  introduce  a 
nonzero  on  the  subdiagonal  of  B.  Similarly,  each  column  transformation 
returns  B  to  triangular  form  while  introducing  a  new  nonzero  element  on 
the  subdi agonal  of  T.  Delaying  the  application  of  the  right  transfor¬ 
mations  to  1  ensures  us  that  the  new  nonzero  element  produced  will  net 
affect  future  row  transformations  on  T,  i.e.,  T  will  remain  in  V.  According 

+  A  f  '  «»•  r\-f*  ssf*  mnl  aw  **%  ft  4-  *•  "?  "  1  S'?*;! 

vv  VltV  .  un  V/A  aocwwxa  VXWIi  VA  lllUXUA^AAVaUAVIi  VA  lUaUA  XI.W.O  UliC  UCxaj  xc  xtgax  • 

In  summa *y,  the  explicit  algorithm  for  each  main  iteration  step  is  given  by: 
Set  T  to  A  -  >B  for  i  =  1,  2,  ...,  n-1. 


1)  Find  B.  to  stably  zero  t.  .  .  and  set  T  to  L.T  and  B  to  L.B. 
i  i+I, i  i  i 
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2)  If  i  >  1,  set  T  to  TMi_1  . 

3)  Find  M.  to  stably  zero  b.  ,  .  and  set  B  to  34. . 

i  1+1,1  l 

Set  T  to  TM  . . 

n-x 

Set  A  to  T  +  \l. 

As  in  the  LR  algorithm,  the  subdiagonal  elements  of  A  should 
become  small  and  approach  zero  at  a  rate  determined  by  the  ratio  of  the 
eigenvalues-  By  using  a  shift  X  ,  we  hope  to  hasten  the  process. 

There  is  only  one  major  drawback  to  the  algorithm  just  described: 
it  is  potentially  unstable.  If  the  shift  X  is  large  relative  to  the 
size  of  the  elements  of  A,  information  needed  to  find  future  small 
eigenvalues  can  be  lost  when  T  is  explicitly  formed.  This  vill  occur 
when  the  shift  is  computed  from  the  lower  2x2  submat rix  of  AB  "  ?.nd 
much  smaller  elements  appear  in  the  bottom  of  B  tfc?n  in  the  top  «:i  2- 
The  following  example  indicates  the  deter  xoraiior.  chat  can 
occur  with  the  explicit  algorithm.  The  relative  residual  is  the  quantity 

Ijg.Ax  -a.Rxjj/(|B.|  !|A|!ro  +  [aj  fj  I 
00 

th 

where  pi  is  the  i  diagonal  element  of  the  fi.ua  1  triangular  mat.  lx 
LBM  and  cr  is  the  i  diagonal  element  of  the  final  vuia.jgular  matrix 
IAM.  The  significance  of  this  quantity  is  that  we  have  really  solved 

X  u 

the  problem  p(A+E)x  =  a(B+F)jc  .  The  i  eigenvalue  is  given  by 
This  problem  was  done  on  an  IBM  3&0  machine  in  double  precision. 


* 
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*->c *5*^5?^»^7'rV2*VW-i.  -* 


"V.-^  v?A  ^CS^J!y^3^| 


Eigenvalue 

- . 699991999943633H*1021 

6 

- . 1142&7 0204249628*10 

0.0 


Relative  Residual 
•7242069F7 6452l8l*10-17 
.8643463693672228*10-2 

.2666666666666667 


According  to  these  results,  zero  is  an  eigenvalue  of  the  problem, 
which  contradicts  the  fact  that  A  is  nonsingular,  and  in  two  instances 
the  relative  residual  is  so  large  that  their  corresponding  eigenvalues 
must  solve  a  problem  which  cannot  be  considered  close  to  the  original 
problem. 

The  instability  mentioned  above  can  be  avoided  if  T  is  never 
formed.  We  will  now  describe  an  implicit  algorithm  which  works  with  A 
and  B  directly.  When  this  new  method  was  applied  to  the  above  example, 
the  following  results  were  obtained: 


Eigenvalue 

Relative  Residual 

- .699992  0003 4 9999*102 1 

.46440675197 8657*10"17 

- .  14  00015 183 15642*107 

.605 69507 054747 8*10"17 

.183673576484812 

-l6 

.149715391676453*10 

The  small  relative  residuals  indicate  that  the  eigenvalues  solve 
a  problem  which  is  close  to  the  original  eigenvalue  problem. 

We  note  that  with  the  standard  eigenvalue  problem  Ax  =  lx, 
Gershgorin’s  theorem  (19)  assures  us  that  computing  the  shift  from  the 
eigenvalues  of  the  lower  2  by  2  of  A  will  not  give  us  a  shift  larger 
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than  the  norm  of  A. 


In  our  description  of  the  new  implicit  algorithm  transformations 
will  be  denoted  by  and  but,  as  we  shall  prove,  the  implicit  and 
explicit  algorithms  are  essentially  equivalent  in  the  absence  of  roundoff 
error,  and  except  in  one  instance,  the  L^'s  of  the  implicit  algorithm  !are 
the  L^'s  of  the  explicit  algorithm.  The  same  is  true  for  the  M^'s'. 

The  implicit  algorithm  begins  by  fonning  the  same  that  was 

*  i  t 

formed  in  the  explicit  algorithm  from  a,,,  b  ,  a„,  and 
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The  matrix  is  applied  to  A  and  B  and  obviously  the  same  nonzero 
element  is  introduced  in  the  (2,1)  position  of  B  as  in  the  explicit 
algorithm.  is  again  formed  so  that  BM^  is  in  T  ,  but 
this  time,  is  also  applied  to  A.  At  this  point  A  and  B  look  like 


A  B 


X  X  X  X 

X 

X  X  X  X  X 

X  X  X  X 

X 

0  X  X  X  X 

X  X  X  X 

X 

0  0  X  X  X 

0  0  X  X 

X 

0  0  0  X  X 

0  0  0  x 

X 

0  0  0  0  X' 

W<*  now  select  from  so  that  L-^V  is  in  V.  When  L,,  is 
applied  to  B,  a  new  nonzero  element  is  introduced  in  the  0,2)  position 
of  B  which  is  then  annihilated  by 

In  general,  row  transformations  return  A  to  upper  Hessenberg 
form  and  introduce  a  nonzero  element  on  the  suodiagonal  of  B. 

Column  transformations  return  B  to  upper  triangular  form  and 
produce  a  nonzero  element  on  the  second  subdiagonal  of  A.  In  contra:  t 
to  what  occurs  in  the  explicit  algorithm,  in  the  implicit  method,  ce1v  n 
transformations  are  applied  simultaneously  to  both  matrices.  T”  mov< 
detail  each  iteration  of  the  implicit  algorithm  is  given  by: 


I 


,  1)  Set  V  to  a1]L  -  X:bi;L  and  6  to  agl  . 

2)  If  |  6 1  >  |  V  |  ,  is  the  element  of  using  pivoting  with 

7]^  = -y/6 ;  otherwise  is  the  element  of  £^  without  pivoting  with 
•  1]^  =  -6/Y.  Set  Alto  L^A  and  B  to  L^B.  Set  i  to  1. 

3)  Find  M^,an  element  Of  7^,  to  stably  zero  b^+1  ^  and  set  A 
to  A>L  and  B  to  BM^ .  If  i  =  n-1,  stop. 

4)  Set  i  to  u+l. ,  Find  L^,an  element1  of  £^,  to  stably 

•zero  a..,  .  ,  and  set  A  to  L.A  and  B  to  L.B.  Return 
i+1, l-L  -l  -l 

to  3* 

s  •  2  2 

For  the  implicit  method,  about  2n  multiplications  and  2n 

l 

additions  are  required  per  iteration.  If  eigenvectors  are  also  requested, 

2  !  2 

n  multiplications  and  n  additions  must  pe  spent  to  accumulate 

2  2 

the  M's.  .In  contrast,  the  QZ  algorithm  requires  additions,  13n 

2 

multiplications  and  3n  squkre  roots  per  iteration,  and  8n  additional 

multiplications  and  additions  if  eigenvectors  are  requested.  However, 

it  should  be  pointed  out  that  to  keep  the  arithmetic  in  the  real  domain 

for  real  matrices,  each  QZ  iteratiort  is  a  double  step.  Thus  a  fairer 

comparison  might  be  to  compare  one  QZ  iteration  tq  two  LZ  iterations  and 

to  keep  in  mind  that  even  for  real  matrices  LZ  uses  complex  arithmetic. 

For  complex  matrices  a  single  shift  version  of  QZ  is  probably  preferable 

to  a  double  shift  .version  of  QZ.  A  single  shift  QZ  iteration  would 

require  6n  multiplications  and  6n  additions  and  2n  square  roots  and  an 
2  2 

extra  3n  multiplications  and  3n  additions  if  eigenvectors  are  requested. 
These  statistics  seem  to  indicate  that  the  LZ  method  is  the  more  efficient 
than  the  QZ  method  for  complex  matrices. 
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It  is  also  interesting  Lo  compare  each  iteration  of  the  LZ 

algorithm  with  each  iteration  oi'  the  standard  LR  algorithm  as  given  in 

2  2 

[  17]  •  The  standard  LR  requires  n  additions  and  n  multiplications 

2  .  2 

per  iteration,  and  n  more  multiplications  and  n  more  additions  if  eigen 
vectors  are  requested.  Thus  the  basic  LZ  algorithm  does  twice  as  much 
work  per  iteration  as  the  LR  method,  but  only  3/2  times  as  much  work 
when  the  accumulation  of  matrices  to  obtain  eigenvectors  is  considered. 

The  operation  counts  given  above  can  be  summarized  as  follows: 


OPERATION  COUNTS  PER  ITERATION 


Without 

Eigenvectors 

Eic 

With 

;envectors 

+ 

X 

Square 

roots 

+ 

X 

Square 

roots 

*LZ 

2n2 

2n 

0 

3n2 

3n2 

0 

Double  QZ 

13  n2 

Hr? 

3n 

2  In2 

2  In2 

3" 

♦Single  QZ 

6n2 

'2 

on 

2n 

9n2 

9n2 

2n 

2 

2 

~  2 

„  2 

*LR 

n 

n 

0 

2n 

2n 

0 

Always  uses  complex  arithmetic. 


The  operation  counts  reported  for  the  double  QZ  algorithm  are 

those  given  in  [10].  If  the  left  eigenvectors  of  the  problem  with 

transposed  A  and  B  are  computed,  as  opposed  to  the  right  eigenvectors 

of  the  original  problem,  then  the  left  hand  transformations  must 

2  2 

be  accumulated  and  only  l8n  additions  and  l8n  multiplications  are 
required  per  iteration. 


In  this  section  we  will  prove  that  the  explicit  and  implicit 
schemes, described  in  the  previous  section, generate  and  use  the 
same  transformations. 


Theorem.  Let  and  M  represent  the  transformation  in  the  implicit 

method  and  L^.  and  represent  the  transformations  in  the  explicit 

method.  If  the  next  iterate  is  A1  and  a'  is  nonzero  for  j  <  i, 

0  >0-1 

then  for  J  <  i,  L  =  L  and  M.  =  M.. 

do  0  0 

I*oof.  The  proof  is  by  induction  on  j.  By  construction  L±  is  equal 


to  L1  and  ^  is  equal  to  M^.  We  assume  that  and  for 


k  <  j,  and  let 


A*  = 

Xj  .  i  •  • 

J-i 

.L-jAM^ 

.  .M .  , 

o-i 

B*  = 

L.  *»  •  • 
0-1 

.  .M.  , 

o-i 

T*  = 

L  •  ■*  *  • 

•  L,TM,  . 

.  .M 

o-i 

1  1 

0-2 

0 


A* 

T* 

B* 

X  X  X  X  X 

X  X  X  X  X 

X  X  X  X  X 

X  X  X  X  X 

X  X  X  X  X 

0  X  X  X  X 

0  X  X  X  X 

0  0  X  X  X 

0  0  X  X  X 

1 1  j  •  >;  >; 

"  "  1  t. 

x 

O  OH  X  K 


U  U  O  X  X 


0  0  0  0  X 


Because  T  =  A  -  XB,  we  must  have 


A*  =  T*M.  ,  +  XB* 
<3 


which  implies  that 


A**  =  L.A*  =  L.T*M.  ,  +  XL.B*. 

3  J  3-1  0 

*  ** 

We  know  that  L.T*M.  .  is  in  V.  Since  L.B  is  also  iny  ,  A 
J  0-1  J 

must  also  be  in  y  .  But  in  the  implicit  method  L  is  the  transformation  from  £ 

which  returns  A*  to  upper  Hessenberg  form  and  if  either  or 

a*  is  nonzero,  there  is  only  one  element  belonging  to  £.  which 

j,J-l  J 

can  accomplish  ttiis.  Since  transformations  to  A  occuring  after  do 

not  affect  the  j-ls^  column  of  A,  the  element  a7.  .  ,  is  nonzero  only 

J  5  J  X 

if  a^5*.  i  is  nonzero.  Since  ^  -t  is  zero  only  if  both 

j,j-a 

a*  and  a*.  ,  .  ,  are  zero,  the  hypothesis  to  our  theorem 

0,0-1  0+1, o-i 

assures  us  that  there  is  only  one  transformation  from  £..  which  couxd 

return  A*  to  an  element  of  y.  Since  both  and  belong  to 

we  know  they  must  be  identical.  By  construction  NL  and  ML  must  also  be 

identical  and  therefore  we  have  proved  our  theorem  by  induction. I 

If  a'  is  zero  then  we  have  no  assurance  that  row  and 

0,0-1 

column  transformations  subsequent  to  in  the  explicit  and  implicit 

algorithms  are  identical,  but  this  is  of  little  consequence.  In  fact, 
the  best  policy  in  both  algorithms  is  that  as  soon  as  a  permanent  zero 
is  detected  on  the  subdiagonal  of  A  ,  then  the  iteration  should  be 
discontinued  and  work  begun  on  a  problem  of  lower  dimension.  If  in 
both  methods,  this  policy  were  adopted,  then  the  algorithms  would  be 


equivalent  up  to  roundoff  error. 


IV.  FINDING  THE  EIGENVECTORS 


After  A  and  B  have  been  reduced  to  triangular  form,  the  eigen¬ 
vectors  can  be  easily  determined.  Let  ay  and  denote  the  diagonal 
elements  of  the  triangular ized  A  and  B  and  let  y  denote  the  corres- 
ponding  eigenvector,  that  is 

(s/  -  “iB)  it  ' 0  • 

The  components  y  of  yA  can  be  obtained  by  solving  this 
triangular  system  as  follows: 


yij  =  °  for  1  <  ** 
y13  *  1  for  i  «  A 


f±i  ~  ^daik  “  ai  bikJ  ykj 

J  i  J  l  k=i+l 


for  i  =  j-l,J-2,...l. 


The  j  eigenvector  of  the  original  system  can  be  found  by 
multiplying  y  by  M. 

If  the  denominator  in  the  above  formula  is  zero,  then  it  is 
replaced  by  macheps  *(||A||  +  ||  B  [j^ ) .  The  denominator  is  zero  when  the 

i  and  J  eigenvalues  are  equal.  If  the  numerator  is  also  zero, 
then  linearly  independent  solutions  will  be  produced.  However,  if  the 
numerator  is  not  zero,  then  y  will  have  large  components  and 
after  normalization  y^  and  ^  will  be  nearly  linearly  dependent. 
This  occurs  when  the  eigenvalue  does  not  have  a  full  set  of  eigenvectors. 


See  Peters  and  Wilkinson[l6]  for  further  discussion 


nvjK  w^T*1 ffwwy^^wt  i  t  w  ^v>  ~- w  w  x?  v^v-'^mTv  >v^  ^  vr«T,^ 
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CHAPTER  TWO  _. 

CONVERGENCE  RESULTS 

In  this  chapter  we  will  prove  several  convergence  theorems,  each 

■  i 

of  which  restricts  the  problem  in  some  way.  A  theorem  might  specify  some 
property  of  the  eigenvalues  themselves  or  characterize  the  matrices  L  and 
M.  Usually  both  types  of  restrictions  are  given..  Many  of  the  theorems 
refer  to  a  nonshifting  or  constant  shifting  version  of  the  I Z  algorithm. 

We  have  found  that  in  most  instances  when  a  shift  policy  has  not  worked, 
the  method  has  been  using  the  same  shift  for  each  iteration.  ,The  chapter  ends 
with  a  partial  listing  of  $  x  5  examples  for  which  the  current  algorithm, 
which  uses  a  solution  of  the  lower  2x2  problem  of  Ax^  =  XBx  as  a  shift, 

i 

will  not  converge  without  the  use  of  an  appropriate  ad-hoc  shift.  < 

We  » ill  use  Parlett's  [  11]  terminology  and  say  thht  a  matrix  is 
an  Unreduced  Hessenberg  matrix  (UHM)  if  it. is  an  upper  Hessenberg  matrix 
and  none  of  its  subdiagonal  elements  is  zerp.  To  simplify  our  proofs 
we  will  assume  we  are  working  with  the  algorithm  in  its  explicit  form. 

The  matrices  will  be  n  x  n,  and  unless  stated  otherwise,  we  will .assume 
that  A  is  a  UHM  . id  B  is  triangular.  For  uniformity  we  will  assume  that 

j 

the  k  iteration  in  the  algorithm  is  given  by 

1)  Find  a  shift  X  .  ! 

K 

2)  Form  Tk  •  Afc  -  ».rBr. 


i)  Find  such  that  L^T^  is  triangular. 

U)  Find  M  such  that  L.RM.  is  upper  triangul 

K  K  K  K 

5)  Set  B,  .  =  L.  3  M,  -  A.  .  =LTR  + 

'  R+l  \  X  K  k-Kl  K.K  K  K  K+l 
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1^  is  a  product  of  matrices  n  n.2 . 1  an<*  s  a 

product  of  matrices  M^.  ^  g . ^knl  where  9ach  ^  is  an  1 

element  of  £^,  and  each  ^  is  in  ^  as  described  in  the  notation  sec 

tion  of  the  introduction.  In  this  chapter  the  multiplier  in  ^  will 
be  denoted  by  T|k  ^  and  the  multiplier  in  ^ :  will  be  denoted  by 

^kji  ' 

In  many  of  our  theorems  we  will  drop  the  iteration  subscript  and 
designate  the  matrices  etc.  by  A,  B,  and  the  matrices  A^+1,  Bj^, 

etc.!  by  A',  B*.  Tiie  matrices  ^  and  ^  will  be  denoted  by  and 
respectively,  and  their  corresponding  multipliers  by  1)^  and  p,  . 
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I.  DEFINITION  OF  CONVERGENCE 


>  • 


[  » 


In  general  the  algorithm  will  be  said  to  be  "convergent"  if,  as  k 
approaches  infinity,  one  of  the  elements  on  the  subdiagonal  of  Ak 
approaches  zero  or  if  for  some  finite  k,  one  of  the  elements  on  the 
subdiagonal  of  is  zero.  Because  of  the  interrelationships  that  exist 
between  A^T  ,1^,  and  ,  we  will  also  regard  the  algorithm  as  convergent 
if  one  of  the  following  conditions  is  satisfied: 

1)  As  k  increases,  one  of  the  elements  on  the  subdiagonal  of 

T  approaches  zero;  or  for  some  finite  k  ,  one  of  the  elements 
k 

on  the  subdiagonal  of  equals  0  . 

2)  For  a  fixed  i,  as  k  increases,  L  .  approaches 

k,i 

the  identity  matrix. 

3)  For  a  fixed  i,  as  k  increases,  ^  approaches 
the  identity  matrix. 

The  reason  for  the  first  criterion  is  that  in  the  explicit  algo¬ 
rithm  the  subdiagonal  elements  of  A^  and  T^  are  identical.  The  reason 

(k} 

for  the  second  criterion  is  that  if  t:  ;  .  is  zero,  then  L.  .  will 

1+1,1  K,1 

he  tli*..*  identity  matrix,  and  tlu:  reason  J'or  the  third  condition  is  that 

if  M,  .  j::  I,  then  will  be  zero. 

k,i  ’  i+l,i 

It  should  be  emphasized  that  when  we  say  IZ  converges,  we 
mean  that  the  problem  can  be  divided  into  two  subproblems  of  lower 
dimension.  We  do  not  necessarily  mean  the  algorithm  can  find  all 


II.  SINGULAR  MATRICES 


o 


If  a  is  zero,  M.  .  is  the  identity  matrix  and  b.  ..  .  ..  is  0  . 

If  b  is  zero,  then  M.  ^  permutes  the  itb  and  i-lSt  columns,  and 

b.  ,  .  .  is  also  0  .  If  b  is  nonzero  and  a  is  nonzero,  then 

l-l, l-l 

bi-l,i-l  iS  either  a-^i_ia/^Tli_ib^b  or  b“(Tli_ib/Tli_1a))b  •  In 
either  case,  it  is  zero  if  the  arithmetic  is  exact.  Again  future  trans¬ 
formations  on  B  during  this  iteration  can  not  affect  the  zero  in  the 
(i-1, i-1)  position. 


We  see  that  with  each  iteration,  a  zero  on  the  diagonal  of  B 
moves  up  one  row.  Within  n-1  iterations  it  must  reach  the  (1,1) 
position,  at  which  point  the  algorithm  must  converge  in  one  iteration.! 

The  following  lemmas  consider  the  case  in  which  B  is  nearly 
singular.  The  quantity  £  is  assumed  to  be  a  small  number  relritive  to 
the  norm  of  B. 


[aim  1:  If  b^  ^  -  £  for  i  >  1  and  L._^  involves  pivoting 

then  l|h,i_Jji_}||  <-  || q I 


Proof:  If  B  represents  the  matrix 


L.  ,l».  ,  .  .  •  •  .M. 

i  - 1  i  •  i  i  i  • 


I  hell  l>  .  .  t  li  n.  Thu:.  l/  will 

i  - 1 ,  i  - 1  i-1, i-1 

he  £  orui_if3  and  sincej]  <  1,  j|  b\  ||  be  £j|g|  . 

1--L  >  1-X 


■P6™.3.  ***  b  .  -  €  -for  i  >  1  and  L.  does  not  involve 

1“ X 


pivoting,  then 


'i  b*i-l,i-lif  * 

i-i 


t 


Proof:  If  B  is  the  matrix 


L.  ,  L.  /-,•  •  •  •L1 BM. .  .  .  •  M . 

l-l  i-2  1  1  1-2 


then  we  may  write 


t)  =  a 

i-l,i-l 


b  i,i-l=  ^i-l3 


b  =  b 

D  i-l,i 


*1,1  =  ^i-lb* 


If  Mi  ^  does  not  involve  pivoting,  then 


i-l,i-l 


« +  ii.ib 

ag _ _ € _ 

€  +  1U-ib  =  Hi-1 


^i-ia 

[rnrjj 


and  hence 


Hb'i-i,i-3!i  -"vT11 


If  ^  does  involve  pivoting  then 


i- !  ,i - 1 


l.  -  ('  •**!  j  _  , l* )•'•/  (llj _  l,‘ ) 


_ 6 _ 

1  i-1 


The  previous  two  lemmas  indicate  that  small  elements  on  the 

r.r  TJ  r.-roo-r>  im  fVio  Hioimnal  rtf*  R  The  npvt  1  Pnwnn  gives  US  an 

V(X  VA  ^  ^  v**>'  “'O - -  -  *  0 

idea  of  what  occurs  when  the  small  element  reaches  the  too  of  B. 
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.00 


Lemma  3 :  If  ^  =  f  and  L^  1  and  M.  .  .  do  not 


\+i,] 


involve  pivoting  for  i  <  j ,  then  there  are  constants  c .  and  d .  such  that 

3  J 


b(kTJ)  =c  .f  end  a'“'*'  =  d.£ 


(k+j)  _  .  j 

2,1  Y 


1,1  J 

for  4  >  0. 

Proof:  Our  proof  is  by  induction  on  j.  For  j  =  1,  the 
hypothesis  to  our  lemma  implies  that 


e  n 


lk,l 


k,l 


~b^  7”n  ^00" 

2,2  ''k,!  1,2 


s  c  6 


which  means  that 


Y  '  '  (1  +  eb£*  >  *  '3.6 

and  if  ^  does  not  involve  pivoting 

a(k+l)  =  _  ,M  .  T1  t(k)  x  _ 

2,1  ,2  '‘k,l  1,2  -  dl^ 

and  if  2  does  inv°lve  Pivoting 


^l>  -  -u!k).  -  d «. 


If  we  assume  our  lemma  Is  true  for  j,  then 


T)  =  -d  (0/  +(xt+j) 
k+j,l  f  '  "11 


which  means  that 


^k+j,l 


d  .c 
3  3 


.HI  (k+j) 

f  '  %1 


"(k+j)  +  p  b(k+j) 


2,2 


k+j ,1  1,2 


=  c€ 


j+1 


33 


so  that 


-tr1  ■  6(cj  +  «&'  > * ^ 

if  L  .  .  does  not  use  pivoting 

K+J  ,2 

(ktj+1)  =  J+l/  t(k+J)  .  -  .  t(k+^ 

a2,l  c€  (  *2,2  +  ^,1  1,2 

-0+1 


s  do+ie 


If  L^.  2  uses  pivoting 

(k+j+l) 

2,1 


-  <***  t'^3)  -  dltl^+1.  | 


0+1 


Our  lemma  did  not  state  the  size  of  the  elements  c  and  d. 

J  J 

and  there  is  no  guarantee  that  they  will  be  small. 


Our  next  theorem  is  a  counterpart  'f  Theorem  i.  We  note  that  B 
might  be  singular. 

Thee  rem _g :  If  T  is  singular,  i.e.  X  is  an  eigenvalue  of  the 
problem  Ax  =  \Sx,  then  the  IZ  al’orifhn  converges  in  one  step  in  exact 

A.' 

arithmetic. 

Proof:  M'Im-  first  n- 1  ••niumns  of  T  jmist  l>o  li hourly  independent 

>»r  else  some  subdiagonal  ••lemon'  of  T  would  be  zei’o  thus  implying  T  is 

not  a  UHM.  The  algorithm  constructs  a  nonsingular  matrix  L  such  that 

IT  =  R,  an  upper  triangular  matrix.  Since  the  first  n-1  columns  of  T 

are  linearly  independent,  the  first  n-1  columns  of  R  must  also  be. 

Similarly,  since  T  is  singular,  R  must  also  be.  This  means  that  the 

last  column  of  R  may  be  written  as  a  linear  combination  of  the  first 

n-1  columns,  and  because  the  last  component  of  each  of  the  first  n-1 

columns  of  R  is  zero,  the  last  component  of  its  n  column  must  also 

T 

be  zero.  Hence  the  last  row  of  R  is  6  ,  the  null  vector. 

3U 


U 


Now  in  the  next  step  of  the  algorithm  we  construct  M  so  that 

B '  =  LBM  is  upper  triangular  and  set  A '  =  T .'  +  kB '  where  T  '  =  BTM. 

Since  the  last  row  of  R  is  0  ,  the  last  row  of  T'  must  also ! he  0  , 

and  hence  a7  ,  must  be  zero.  I  .  i 

n;n-l 

Our  theorems  and  lemmas  indicate  that  in  the  future  we  can  safely 
ignore  problems  where  either  A  or  B  is  singular.  However,  we  would 
like  to  include  singular  cases  in  the  next  few  theorems  wherever 
possible^ because  many  of  these  theorems  not  only  guarantee  convergence, 
they  also  give  some  hints  about  rates  of  convergence. 
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III.  Cl«01iAi.  COIWKKGKNCK 
.  ~  -  -  •; — 


In  uhis  section,  we  phall  prove  several  global  convergence 
theorems.  Most  of  our  results  refer  to  a  constant  shifting  algorithm,  although 
in  practice  our  shifts  are  usually  different  for  each  iteration.  How¬ 
ever,  we  have  often  found  chat  when'  a  shift  policy  has  not  worked,  the 
same  shift  is  being  used  for  each  iteration.  Thus,  our  theorems  do 
have  practical  significance  although  they  do  not  refer  to  any  actual 

'  i 

implementation  of  I Z.  ;  : 

We  begin  by  proving  four  similar  theorems.  They  all  consider 
the  sequences  of  matrices  M  -  W  *  W  -  h\  ’  M  “d  M 

l 

where  ^  =  \  f°r  some  scalar 


\+i  VHA 


B, 


k+1  "  I‘kBkMk 


sk  -  Vk 


=  A 


Bl  =  B 


and  1^  and  are  nonsingular  for  all'  k  and  Sk  and  are  upper 

.  '  i 

triangular  for  all  k. 


Tli'  or  m  . :  ll  B  is  non.. angular. 


Vk=Lx 


V1* 


and 


then 


U,  =  S,  B,"1  S.  B*1,  . 
k  k  k  tc-i  k-i 


•SlBi 


V,U.  =  (A  -X.B  )B_i  (A  -X,  .  B)  B_1 - (A  -X.B  )B_1 

rC  K  K  '  K“l  X 


j5o 


Vk  ul"l  %  =  ((A  -  PB)B_1)^e^*  I 


Theorem  4.  If  A  -  \.B  is  nonsingular  for  all  i,  then 

Uk“\_1  =  B  (A  -  ^B  )_1B  (A  -  XgB  )_1  ....  B  (A  -  XfcB  )_1. 


The  proof  of  the  above  theorem  parallels  that  of  the  previous  theorem. 


Corollary  4.1.  If  the  conditions  of  Tneorem  4  are  satis- 

-T 

fied  and  \,  =  p  for  all  k,  then  V.  e  is  in  the  direction  of 
k  5  k  ~n 

((A  -  pB  )  B  )  en. 

Proof.  By  theorem  4 


U'V  =  [B(A  -  pB)_i)k 


which  means  that 


\T\T  =  ((AT  -  pBT)'1BT)k  . 


Since  is  the  product  of  upper  triangular  matrices,  the  matrix  P^sU^  must  be 

(k) 

lower  triangular  which  implies  that  P,  e  =  p  '  e  and  hence 
^  c  k~n  rn,n  ~n 

p«  -  «*T  -  »BT)-Wen.  I 


If  the  matrices  described  above  represented  the  matrices  in  a 
version  of  the  LZ  algorithm,  which  did  not  allow  row  interchanges,  then 
for  all  k,  the  matrix  would  be  unit  lower  triangular. 


t 


G 


Theorem  5.  If  b  is  nonsingular. 


\  -  E5i-i!«1- 


•B-1  Sx 


wk  -  . \ 


Vk  =  bi'3<ai  -  xkBi)Bi1(Ai  -  Xk-1B1}-  •  •  •  Bi>i  ■  xiBi> 


Proof.  Since 


Tk  -  *k-i(*k-i  -  xkBk-i)Mk-i 


“k"1  ■  C A-1C1  > 


we  have 


°  -  ‘AAi 


-  wk-iBi1<Ai  -  xkBi>wk-i 


which  implies 


VA\  *  BI1<A  -  xkB  )wk-i 


“A "  '►-iWlAVi 


-  “k-A^V^Vk-i 


TT  n"lm  TT 

"k-A  Aku: 


lm  TT 

;  xkuk-l 


“B  <A  -  A-A-l 


=  B-1  (A  -  XkB  )ET1(A  -  Xk_1B  ). 


.  -  >  .3 


Corollary  5.1,  If  the  hypothesis  of  Theorem  5  is  satisfied  and 

T  k 

\j.  =  p  for  all  k,  then  is  in  the  direction  of  (B'^A  -  p3))  e^. 


Proof.  By  Theorem  5 


WkUk  =  (B-1(A  -  pB))  . 


Since  Uk  is  a  product  of  upper  triangular  matrices,  Uke1  = 


u«  e 
1,1  ~1 


which  means  that 


Wk  *i*i  fl  =  (B’1(A  '  PB»\-  1 


Theorem  6.  If  A  -  \kB  is  nonsingular  for  all  i  ,  then 


U^W'1  =  (A  -  11B)"XB(A  -  X2B)~XB.  .  .  .(A  -  XkB)"xB 


The  proof  of  this  theorem  parallels  that  of  the  previous  theorem.  Note 

i  1 j_  4.1.  *  „  Ja/%'*  v\/n4-  AAAtimA  f  V««+  D  •?  o  nAnei  nmil  or 

oiidu  unxc  oucvi  mi  us/Ou  uwv  aoovunc  • 
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Corollary  6.1.  If  the  hypothesis  of  Theorem  6  is. satisfied 

and  X,  =  p  for  all  k  then  W.^0  is  in  the  direction  of  (BT(A^  - ' p )  ^)ke 
x  k  ~n  •  *  v  ~ n 


Proof.  By  Theorem  6 


which  means  that 


UvX1  =  ((A‘  -  pB)"^)* 


wv\T  =  (3T(A  *  pB)'T)k  . 


Since  is  a  oroduct  of  upper  triangular  matrices,  Q^sU~  is  lower 

(k)  ! 

triangular  and  en  equals  q^‘  £  en-  This  in  turn  implies  that 

-T  (k)  ,DT,AT  T^v-lxk  ® 

W,  qw  e  *=  (B  (A  -  pB  )  )  e  . 

k  ^rijn  „n  .  K  .  '  '  ,„n 


If  we  again  relate  the  above  theorems  to  the  LZ  algorithm,  then 

would  be  unit  lower  triangular  if  column  interchanges  were  not  permitted. 


If  D  is  a  diagonal  matrix  with  diagonal  elements  d1,dp,...dn 

then  the  Moore-Penrose  pseudo  inverse :of  B,  is  the  diagonal  matrix , denoted  by 

+  * 

D  ,  v.ith  elements  z.,z0,....z  where 

1  2  n 

z i  =  l/d^  if  d^^  i$  nonzero’ 
and  z^  =  0  if  is  zero. 


Ve  will  denote  by  i  the  projection  matrix  U+U. 


i  for  i  =  j  and  d^  is  nonzero. 


0  elsewhere 


WMMl  ff!M  W  tut-  ypnxyj rejTOfygw 


r 


JTO!®S8^^?^®^R^9S^^  ^S^^S^S&RSfWWSEWSWI  STi^ir  v.'Vs^vP  *~>m&1l  -  v>^>—  u 


r% 


\  The  next  main  theorem  is  a  modification  of  one  that  appears  in 
Wilkinson  [ig]  and  Parlett  [if]  among  other  places.  Our  version  of  the 
theorem  extends  their'  results  to  cover  singular  matrices .  If  a  'matrix 
X  is  said  to  have  an  LU  decomposition,  then  there  exists  an  upper 

*  I 

triangular  matrix  U  fend  a  lower  triangular  matrix  L  such  that  X  =  LU. 
This  decomposition  is  unique. 

Theorem  7 .  If  F  is  a  matrix  with  (eigenvalues  of  distinct 
modulus  satisfying 

1  •  ; 

Hii>  |d2|>  .  .  .  .  >  |dj>o, 

and  ii  1'  can  bo  written  as  XI)  Y  where  .  Y  =  X  j  J)  -  is  the  diagonal 
matrix  of  eigenvalues  and  Y  has  an  LU  decomposition  and  x  has  an 
Kj  decomposition  L  U  ,  then  the  lower  triangular  factor  of  the  LU 

A  A 

•  S 

decomposition  of  Fk  goes  to  L  $  as  k  -♦  ® . 


:| 

I 

r 

3 


* 

t 

5 


To  facilitate  the  proof  of  Theorem  7 ,  we  first  present  a 


few  lemlnas. 


I^nmaji.  If  L  is  a  lower  triangular  matrix  and  D  is  a  diagonal 

I 

matrix  whose  elements  satisfy 

'Kl  >  l«al>  •  •  •  >  |dj>  |d  |  =...  =|,y 


Un‘ii 


BsawBVsvrr  ;-w«ksb-ju.  _ 

_  •'  '  ~  - 


ana 


Di.  =  DU)  I). 


Proor.  II-  s  =  Dl„  then  a,  ,  -  0  for  j  >  i  or  i  >  ra 

3  yd 

and  s.  .  =  d.-t.  .  elsewhere. 

•*•>«?  -*•  i  }  J 

Now  the  last  n-ra  columns  of  S  are  2ero  because  for  j  >  m, 


s.  .  =0  for  i  >  m 

1  >J 


s . 


.  .  =0  for  i  <  j,  i.e.  for  i  <  m. 

1  >d  — 


Multiplying  S  by  D  D  zeroes  the  last  n-m  columns  of  S  and  leaves  the 
first  m  columns  untouched.  Since  the  last  n-m  columns  of  S  are 
already  zero  we  have  S  =  SD+D  or 


DL  =  DLD  D. 


Lemma  5 •  Let  L  be  a  unit  lower  triangular  matrix,  and  let 


D  be  the  matrix  oi  Lemma  4,  and 


G  -  DKLD+K; 
k 


then 


G.  -  $  +  K  where  E  ••>  O  as  k  ->rt~ 

K  K  K 


Proof.  If  g is  the  (i,j)  element  of  G,  ,  then 
ifj  K 


■S’-  { 


0  for  i  <  j  or  j  >  m 

L  (d./d.)1^  elsewhere  * 
i»j  1  J 


Tlio  fact  Uiat  |<h/d^.|  <-  J.  for  i  >  j  and  ,j  <  m 

and  ii./d.  -  L  for  i  <  in,  togi:th<-r  with  the  fact  that  L  is  unit  lower 

i  i  — 


'k) 

triangular  implies  that  g^.  =1  for  j  <  m  and 

<]  J 
(k) 

g. /-»  0  for  j  <  i  and  j  <  m  as  k  approaches  infinity.  Thus 

1J  ■ 

Gk  =  $  +  Ej^  where  ->  0  as  k  approaches  infinity.  ■ 

Lamma  6.  If  U  is  an  upper  triangular  matrix  and  D  is 

the  diagonal  matrix  of  Lemma  k -  then 

USD  =  SUD  . 


and 


Proof:  If  G  =  TJ5  ,  then 

g  =0  for  j  >  m  or  i>j 

i»-j 

e  =  u.  .  elsewhere. 


if  H  =  GD,  then 


and 


0  for  j  >  m  or  i  >  j 

u.  .d.  elsewhere. 

i,J  0 


The  formulae  just  given  for  h.  .  are  exactly  those  which 

would  be  given  for  e.  .  where  E  =  UD.  Thus  H  =  E.  Moreover,  the  last 

1,0 

n-m  rows  of  E  are  zero  so  that  i  E  =  E .  Thus 

$E  =  H  or  $TJD  =  USD.  ® 


Wc  can  now  give  the  proof  of  theorem  7 • 

Proof  of  theorem  7 :  Assume  the  eigenvalues  of  the  matrix 
F  are  of  distinct  moduli  and  assume  F  can  be  written  as  F  =  XDX 
where  D  is  the  diagonal  matrix  of  eigenvalues  and  both  X  and  X  1 
have  LU  decompositions.  Let  L  U  be  the  US  decomposition  of 
X  and  let  L  U  be  the  LU  decomposition  of  X  ^  It  can  be  easily  shown 

y  y 

uu 


=  XD  L  D 

y 


+K 


If  G,  =  DKL  D+K,  then  by  Lemma  5  G. 
k  y  k 

k  goes  to  infinity. 


by  Lemma  4 . 


$  +  E^  where  goes  to  0  as 


F*  =  X($  +  E^D*^ 

=  L  U  ($  +  E, 

xx  k  y 

=  L  (#  +  U  E.U~1)U  D^U  by  Lemma  6- 
x  x  K  x  x  y 


Since  ^E^U^1  goes  to  0  as  k  goes  to  infinity,  $  +  U  EjJJ^  goes 
to  f  as  k  ->oo.  Because  $  U^D^U  is  upper  triangular,  the  lower 
triangular  factor  of  approaches  L  i  as  k  becomes  large.  | 


Theorem  8  :  If 

(1)  B  is  nonsingular 

(2)  LZ  is  used  with  a  constant  shift  p. 

(3)  The  quantities  X^-p  for  i  =l,2,...n  have  distinct  moduli 

(4)  Either  no  row  pivoting  is  required  and  there  exists  a  matrix 

X  such  that  (A  -  pB)B  ^  =  XDX  X  where  D  is  diagonal  and 

.1 

both  X  and  X  *  have  LU  decompositions  or  no  column 
pivoting  is  required  and  there  exists  a  matrix  X  such  that 
B  ^(A  -  pB)  =  XDX  ^  where  D  is  diagonal  and  both  X  and  X 
have  LU  decompositions , 


then  LZ  converges. 


Proof:  We  let  F  =  (A  -  pB)B  1  and  apply  Theorem  7  and  find  that 
as  k  increases,  the  lower  triangular  factor  of  the  UJ  decomposition 
of  F*c  approaches  Lx$  .  By  Theorem  3  this  lower  triangular  factor  is 
given  by  Vfc.  If  L$  and  Vk+1  -»  L*  and  Vk+1  =  VkI^  ,  then  L^1 
must  be  approaching  the  identity  matrix  which  means  that  LZ  is  convergent. 

If  F  =  B^(a-pb)  then  by  invoking  Theorem  7  and  Theorem  5 
we  see  that  ^  must  be  approaching  the  identity  matrix  as  k  -»  ®  which 
means  that  12  is  convergent.  1 

The  condition  that  both  X  and  X  1  have  KJ  decompositions 

-1  -1 

is  partially  satisfied  since  both  (A  -  pB)B  and  B  (A  -  pB)  are 

unreduced  Hessenberg  matrices.  Parlett  [14]  has  proved  that  if  F  is  a  UHM, 

-1 

then  there  exists  a  matrix  X  such  that  XJX  is  the  Jordan  canonical 
form  of  F  and  X  ^  has  an  LU  decomposition. 

The  condition  on  the  distinctness  of  the  moduli  of 
the  eigenvalues  can  be  relaxed  somewhat.  Wilkinson’s  [  19]  treatment 
ol  multiple  eigenvalues  for  the  LR  algorithm  can  be  applied  directly 
to  the  LZ  algorithm. 
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IV.  ITERATION  AND  LZ. 


Theorem  8  indicates  that  with  a  constant  shift 
the  LZ  algorithm  converges  for  problems  with  shifted  eigenvalues  of  distinct 
moduli  if  row  pivoting  is  not  required  or  if  column,  pivoting  is  not 

required.  Seme  of  these  restrictions  can  be  weakened  by  investigating  the 
relationships  between  LZ  and  various  iterative  procedures.  This  approach 
has  been  taken  by  Poole  [  13] 5  Parlett  and  Kahan  [  12] ,  and  Buurema  [  1  ] 
in  studying  the  QR  algorithm. 

Lot  us  consider  the  LZ  algorithm  with  constant  shift  p  for 
the  problem  having  eigenvalues  \  ^  and  let  us  denote  \ _  -  p  by  d^ 
and  assume  that  the  fx^j's  are  ordered  so  that 

idj  >  i<u  >••••>  iani  • 


Theorem  9 .  If  any  of  the  following  4  conditions  is  satisfied 
then  the  algorithm  converges. 

1)  | dj  ^  ! a  | ,  e^  is  not  orthogonal  to  the  right  eigen¬ 

vector  of  (A  -  pB)B  1  associated  with  d  and  row  pivoting  is 
never  .lone  to  zero  for  all  k. 

2)  j  d^j  7^  |  <jn  ^5  ,  en  is  not  orthogonal  to  the  left  eigen¬ 
vector  of  B(A  -  pB)  \  associated  with  V^n  and  row  pivot 'ng 


is  never  done  to  zero 


t^  ,  for  all  k. 
n,n-l 


3)  |  d  J  t  i  r‘^1  ,  e^  is  not  orthogonal  to  the  right  eigen¬ 
vector  of  B  1(A  -  p B)  associated  with  and  column  pivoting 

(k) 

is  never  done  in  zeroing  b^  '  for  all  k. 
h)  j d  j  ?  j  1| ,  is  not  orthogonal  to  the  left  eigen¬ 

vector  of  (A  -  p  B)  ^B  associated  with  l/dn  and  column 


pivoting  is  never  done  to  zero  b 


n,n-l 


for  all  k. 


Proof.  The  proof  of  (l)  entails  looking  at  the  power  method 
for  solving  (a  -  PB)x  =d  Bx  given  by 

1)  Set  Xq  to  e^  and  k  to  0. 

2)  Find  y  such  that  By  =  x^  . 

3)  Set  z  to  (A  -  pl)y  . 

*0  Set  x^  =  z/||z||  . 

If  -  x^ll  *s  small  then  stop; 

otherwise,  set  k  to  k+1  and  return  to  2. 

If  |d^|>|d2|  and  e^  is  not  orthogonal  to  the  lei t  eigenvector  of 
(A  -  pB)B  ^  corresponding  to  d^>  then  the  power  'r.ethod  converges  and 
x^  will  be  that  eigenvector.  We  note  that  x^  is  in  the  direction 
of  ((A  -  pB)B  )  e, ,  which,  according  to  Corollary  3.1,  is  also  the 
direction  of  V^e^.  Thus  if  the  power  method  converges,  it  should  be 
clear  that 


V.  ,  e .  =  c  V  e.  +  f 
k+l~l  k  k-1  ~k 


where  c^  is  some  scalar  and  f^  -»  0  as  k  -»  a. 


(2  -  1) 


Since  Vk+1  =  V^1 


and  L 


-1  _  r~l  T-1 


X  •  W 


\,n-l  »here  \,i  ls  in  fi 


which  means  that  IT1  either  has  the  form 


r  1 

9 

0 

X 

0 

^  0 

i 

US 


4 

;  ‘A 

| 

% 


a 

;  I 


$ 


t 


c 


* 


then  either 


Vk+lfl 

Vkfl 

\,lVkf2 

'  (2  -2) 

or 

V  e 
k+lfl 

Vkf2 

\,lVkfl  * 

If  pivoting 

is  never 

(k) 

done  to  zero  tj>  ; 

2,J- 

then  equation 

(2-2)  holds  for  all  k  and  V,  e  is  linearly  independent  of  the  other 
columns  of  Vk,  and  if  equation  (2-1)  is  also  satisfied,  T)k  ^  must 
be  approaching  0  and  ^  must  be  approaching  the  identity  matrix, 
i.e.  the  LZ  algorithm  must  be  converging.  If  equation  (2-1)  holds  but 
the  first  column  of  V.  is  approaching  a  multiple  of  the  second  column 

K 


t 


9 


9 


of  V,  j  then  we  cannot  assert  that  LZ  will  converge.  It  is  possible  ■ 

K 

t 

for  equation(2-l)  to  hold  without  ^  converging  to  the  identity  matrix. 

The  proof  of  (2)  entails  looking  at  an  inverse  iteration  scheme  . 
for  finding  x  such  that  (AT  -  PBT)x  =dBTx  .  The  iteration  scheme  can 
be  summarized  as  follows: 

1)  Set  Xq  to  e  ,  K  to  0. 

m  »p.  T 

2)  Find  z  such  that  (A"  -  pB  )z  -  B  x^  . 

3)  If  z  is  large,  stop. 

k)  Set  x.  to  z/||z|;  ,  k  to  k+1,  go  to  2. 


If  |  dn|  <  jcM^J  and  is  not  orthogonal  to  the  left  eigenvector  of 

B(A  -  pB)  1  corresponding  to  V^n  >  then  the  inverse  iteration  scheme 

will  converge  and  z  will  be  that  eigenvector.  The  vector  x^  is  always 

in  the  direction  (A  -  pB  )  Be  ,  which,  according  to  Corollary  4.1, 

~n 

-T 

is  also  the  direction  of  V,  e  .  Hence  if  the  inverse  iteration  scheme 

ic 

converges,  then  for  large  k  we  may  write 


-T  -T 

Vj\e  =  c,V,  e  +  f, 
k+l~n  k  k  ~n  ~k 


(2  -3  ) 


where  c^  is  some  scalar  and  approaches  0  as  k  approaches  infinity. 


Since  V& 


and  l£ 


-T  "T 

V? 


,T  rT 


\,A,2 . \,n-l  where  V 


is  in  X, 


(C) 


where  X  is  a  dense  matrix 


then  either 


V“T  e  ?  v'Te  + 

vk+l  _n  k  ~n 


1  I  ciV3i  (2“4^ 

lk,n-lj^1  J  *-3 


-T 

V,  T,e 
k+l~n 


I  cATf }  ;  +  \,n-lVkT!n 

(l  ^ 


where  one  c.'s  are  products  of  the  multipliers  involved  in  the  iteration 
J 


s*ep. 


If  pivoting  is  never  done  to  zero  tj^_^  ,  then  equation  (2-4) 

holds  for  all  k  and  the  last  column  of  V^T  must  be  linearly  independent 

of  the  other  columns  of  V^T  .  In  this  case  if  equation  (2-3)  is  satisfied, 

then  Tl,  ,  must  be  approaching  zero  and  LZ  must  be  converging.  If 
k,n-l 
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equation (2 -3 )  is  satisfied  but  the  last  column  of  V  is  approaching 

It 

a  linear  combination  of  the  other  columns  of  V~T  as  fast  as  it  is  converging 
to  an  eigenvector  then  we  cannot  assert  that  LZ  will  converge. 


The  proofs  of  (3)  and  (4)  closely  resemble  those  of  (l)  and 

(2).  For  (3)  we  consider  the  inverse  iteration  method  for  finding  x 

such  that  Bx  =  d(A  -  pB)x  and  use  corollary  5. 1 .  For  (4)  we  consider 

ffl  T  T  -1 

the  power  method  for  determining  x  such  that  B"(A  -  pB  )  x=a.x  and 
use  corollary  6.1 .  | 

Theorem  9  requires  more  elaboration.  It  would  seem  that  because 
the  matrices  and  have  determinant  1  for  all  k,  the  linear  depen¬ 
dence  of  the  columns  of  these  matrices  and  their  inverses  should  not  be 
an  issue.  However,  if  we  look  at  the  singular  values  of  these  matrices 
we  can  get  a  different  picture  of  the  situation.  It  is  possible  for  these 

matrices  to  become  singular  at  the  same  rate  as  one  of  their  columns  is 

-T 

nrmroaching  an  eigenvector.  If  the  last  two  columns  of  V  &i’e  approaching 

k 

the  same  vector,  then  we  cannot  ar.suj.ie  that  pivoting  will  stop  and  the 
algorittim  will  convcrgi .  This  fact  is  brought  out  by  the  following 
2  by  2  example. 


Let  B  be  the  identity  matrix  and  A  the  real  matrix 


a  D 


b  0 


where  | b|  >  | a|  >  0.  If  LZ  is  applied  to  the  above  problem  with 
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a,2/b22  as  the  Policy,  the  shift  is  always  0,  and  A  and  B 

are  left  unchanged  by  each  iteration  and  the  algorithm  does  not  converge. 

-T 

The  least  singular  value  of  is  given  by  the  smallest  eigenvalue 
value  of 


which  is 


where 


t 


-T 

Since  0  <  ja/bj  <1,  |t  |  <  1,  which  means  that  as  k  increases,  V^. 
approaches  a  singular  matiix  and  hence  its  columns  become  multiples  of 
each  other. 

It  is  ironic  that  we  seem  to  require  that  e^  be  nonorthogonal  to 
a  right  eigenvector  of  a  UHM  and  that  e^  be  nonorthogonal  to  a  left 

eigenvector  of  a  UHM.  If  the  situation  were  reversed,  there  would  be 
no  problem,  for  if  en  were  perpendicular  to  a  right  eigenvector  of  a  UHM 
then  the  eigenvector  must  be  9. 

The  last  interesting  fact  about  Theorem  9  is  that  it  suggests 

pivoting  schemes  which  should  converge  although  perhaps  very  slowly. 

If  we  find  that  pivoting  is  necessary  to  stably  zero  t  or  b 

2,1  2 .1, 

we  could  change  the  shift  so  that  pivoting  is  no  longer  necessary.  We 
could  then  view  the  matrices  as  representing  a  new  problem,  and  hope¬ 
fully  the  two  new  largest  eigenvalues  do  not  have  the  same  modulus. 


V. 


COUhTEREXAMPIES 


9 


9 


S 


t 


t 


t 


This  section  contains  a  partial  listing  of  3  X  3  examples  for 
which  tlu-  TZ  algorithm  will  not  converge  when  the  shift  is  the  eigen¬ 
value  of  the  lowest  2x2  subproblem  of  Ax  =  XBx  closest  to 

rss  rv 

a  /b  .  The  counterexamples  share  several  characteristics: 

nn  nn 

1)  The  subdiagonal  elements  of  A  and  the  matrix  B  repeat 

themselves  every  third  iteration. 

2)  Row  and  column  pivoting  are  necessary  at  each  stage  of  the 

process. 

3)  The  shift  is  the  same  for  each  iteration. 

The  first  property  guarantees  nonconvergence.  The  second 
property  is  a  necessary  condition  for  the  first  property:  if  at  some 
st.age  pivoting  is  not  required  to  maintain  numerical  stability,  then 
cycling  will  not  continue  as  before.  The  third  property  indicates 
that  the  constant  shifting  hypothesis  of  Theorem  8  is  realistic  in  terms 
of  actual  computation.  In  fact,  a  constant  shift  may  be  useful  in  ■ 
practice  ass  warning  ce  npnconvergence.  It  is  significant  that  no 
condition  is  specified  for  the  eigenvalues  of  a  problem.  There  are 
counterexamples  in  which  all  the  eigenvalues  are  of  distinct  modulus. 

The  first  class  of  examples  is  the  basic  one.  In  this  case  the 
matrices  initially  look  like 


where  the  following  conditions  hold: 


|a|  <  |h|  jas/gJ  <  |f| 

|as/h|  <  |g]  |aVg{  <  |f| 

[b  1  <  |fj  |bmh/(gf) j  <  jcj 

|Wf|  <  |gj  |bmh/(cf)|  <  \z\ 


lanvg  I  <  M 

|a  |  <  jcj 
|bn^s|  <  { c  | 

lb  !  <  lcl. 


If  B  is  the  identity  matrix  and  A  is  the  matrix 


then  the  above  conditions  are  satisfied. 

For  problems  in  class  1  the  shift  is  zero  for  each  iteration. 
After  one  iteration  the  matrices  are 


After  two  iterations  the  matrices  look  like 


After  three  iterations  the  matrices  return  their  original  forms. 


The  second  class  of  examples  consists  of  those  problems  whioh 


fall  into  the  first  class  after  one  iteration  but  initially  look  like 


* 


* 


s 


* 


t 


t 


t 


* 


ft 


t 


t 


where 

M  <  lci  •  !.y  -  f  dl  <  1hl 

Ip  -  f  g|  <  N  |  (y  -  |  d)s/h|  <  jm|. 

The  first  shift  is  0  and  after  .one  iteration  the  matrices  of  this 
class  have  the  form 


A  B 


where  a  and  b  are  complicated  expressions. 

i 

l 

The  third  class  of  counterexamples  includes  problems  in  which  A 

and  B  are  initially  given  by  1 

A  B  ‘ 


The  shift  is  d/m  and  after  shifting,  ^ 2  anc*  ^3  are  zero •  If- 
this  new  .  problem  is  in  class  2,  then  in  one  iteration  before  the  shift 
has  been  added  back  A  and  B  would  look  like 
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! 


where,  x  and  y  are  again  nbnsiraple  expressions.  Shifting  back  sets 
a,,,,  ,  and  avy  to  d.  From  here  it  should  be  clear  that  the  matrices 
will  repeat  themselves  every. third  iteration  if  their  elements  satisfy 
the  appropriate  magnitude  relationships. 

If  in  the  previous  example,  a.  •  were  initially  zero,  the 

*  ! 

first  shift  would  again  be  d/m  .  Performing  one  iteration  and 
shifting  back  set  the  element  to  d  and  revert  us  to  the  previous 

exarple. 


These  counterexamples  indicate  that  assuming  the  algorithm 
uses  a  constant  shift  is  reasonable  and  that  given 
this  assumption,  the  structure  of  L  and  M  is  important.  The  above 
examples  were  all  constructed  by  assuming  that  pivoting  was  necessary 
to  maintain  stability  at  every  step  of  the  algorithm.  Indeed  if 
pivoting  ceases  at  seme  stage,  then  "cycling",  will  not  continue  as 
before.  It  is  doubtful  that  without  an  analysis  of  a  given  shift 
strategy  we  can  weaken  significantly  the  criteria  given  earlier  for  global 
convergence.  Moreover,  it  is  extremely  doubtful  that  such  an  analysis 
could  give  us  more  than  assymptotic  convergence  results. 
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APPENDIX 


NUMERICAL  RESULTS  AND  FORTRAN  PROGRAM 

The  algorithm  described  in  Chapter  1  has  been  implemented 
in  a  Fortran  program.  The  program  is  designed  to  find  the  eigensystem  for 
complex  matrices,  and  consists  of  two  subroutines  which  must  be  called 
separately.  The  first  subroutine, GELHES,  reduces  A  to  upper  Hessenberg 
form  and  B  to  upper  triangular  form.  If  A  and  B  are  already  in  these 
forms  and  no  eigenvectors  are  required,  then  calling  GELHES  is  unnecessary. 
The  second  subroutine  GLR  finds  the  eigenvalues  and,  if  requested,  the 
eigenvectors  of  the  system.  The  parameters  involved  in  the  subroutine 
calls  are  described  in  the  comments  at  the  beginning  of  each  subroutine. 
Both  subroutines  use  the  subroutine  RABS  to  compute  the  norm  of  a 
complex  number. 

It  should  be  emphasized  that  the  variable  MACHEPS  is  machine 
1-t 

dependent.  It  is  p  where  the  machine  gives  t  digits  in  base  p.  It 
is  set  for  the  IBM  360  double  precision  mode. 

Our  program  has  been  finding  eigenvalues  which  correspond  to 
problems  close  to  the  given  problems:  our  relative  residuals  have  always 
oiv'n  close  to  the  precision  of  the  machine. 

For  our  test  examples  the  total  number  of  iterations  has  been 
roughly  times  the  order  of  the  matrices.  In  general,  for  a  matrix 
of  order  n,  the  time  required  on  the  IBM  $60  model  67  in  Fortran  H,opt=S, 
has  been  about  .75rr^milliseconds  if  eigenvectors  are  computed  and  .hr? 


The  example  given  below  were  generated  using  integer  arithmetic 
by  multiplying  two  bidiagonal  matrices  by  random  nonsingular  transformations.  The 
problems  were  run  on  che  IBt-?  $60  Fortran  G  compiler.  The  relative 
residual  is  the  quantity 

|  -Ax^  ~  I  as 

W  l!AH  +  Kl  ! i Bl L 

th 

where  cy^  and  are  the  i  diagonal  elements  of  the  triangularized 

th 

A  and  B  and  x^  is  the  i  eigenvector. 

In  the  first  example  A  has  rank  4  and  B  has  rank  3  and  their 
null  spaces  intersect.  Since  the  rank  of  B  is  less  than  the  rank  of  A, 
there  is  an  eigenvalue  which  might  be  regarded  as  infinite  because  a 
small  perturbation  in  B  would  yield  a  large  eigenvalue.  Indeed  an 
eigenvalue  of  lO1^  was  found.  The  problem  is  also  "ill  -disposed", 

because  for  any  vector  x  in  the  intersection  of  the  null  spaces, 
any  scalar  \  will  satisfy  Ax  _  XBx  ,  and  may  be  considered 
an  eigenvalue  ol‘  the  problem.  The  example  also  has  three 

genuine, finite  eigenvalues  which  the  algoritlim  war,  able  to  find  accurately 
up  to  the  precision  of  the  machine  despite  the  presence  of  the  two 
spurious  eigenvalues. 

In  the  second  example  there  are  two  double  roots.  The  first 
corresponds  to  a  quadratic  elementary  divisor  and,  as  expected,  is 
accurate  only  up  to  the  square  root  of  the  machine  precision.  The 
second  corresponds  to  two  linear  elementary  divisors  and  is  accurate 
almost  up  to  the  machine  precision. 


Example  1 


rHE  MATRIX  A: 


575.+  -110.1 

205. + 

-460.  1 

-30. + 

456.  1 

40.  + 

-352.  1 

-1G5.+ 

-332.1 

355.+  -356.1 

3  05.+ 

-324.  1 

-070. + 

184.  1 

730. + 

-200.  1 

235.+ 

100.  I 

345.+  -193.1 

500. + 

-440.  1 

-400.  + 

352.  1 

625.+ 

-412.  1 

-100. + 

-16.  1 

210. +  -200.  I 

-215. + 

29.8.  1 

109. + 

-292. 1 

-115. + 

258.  1 

115. + 

24.  1 

-545.+  150.1 

-115. + 

282.  1 

-30.  + 

-276. 1 

-40. + 

222.  1 

115. + 

248.  1 

THE  MATRIX  3: 

GO . +  -387.  I 

33.+ 

504.  1 

-68.+ 

-012. 1 

143. + 

567.  1 

36.+ 

463.  1 

241. +  0.1 

48.+ 

103.1 

28.+ 

3  6.  1 

139. + 

63.  1 

-76 .  + 

-468. 1 

219. +  -533.1 

70. + 

558.1 

-8.+ 

-504. 1 

1G7.+ 

567.  1 

-24.+ 

-72.  1 

46.+  414.1 

96.+ 

-675.1 

-56.+ 

630.  1 

130. ♦ 

-693. 1 

-52.+ 

18.1 

-91. +  189.1 

-112. + 

-207.1 

7G.+ 

306.  1 

-177. + 

-252. 1 

-8.+ 

-378. 1 

TRUE  EIGENVALUE 


I I'F I '!  ITE  EIGENVALUE 


-o.snnnoooooooonooD 


5  -n. 5000000000000000 

4  n.inr.rififiGr»fiGRr»RR7D 

5  ANY  SCALAR 
.\LP«A 


00  + 
00+ 
01  + 


0.5900000000000000  00  I 

■o. sooooononnoooooD  oo  i 
■O.GGC00R66G66G6G7O  00  I 


1 

2 

-* 

j 

4 


-0.1832106495233700 
-0.1300211517428030 
-P.G4237347G326G1SD 
0.1430000000000000 
-0. 14 G3P. 82010523830- 


03  + 
03  + 
03  + 
03  + 
12  + 


•0.20621265635B070D  03 
•0. 1137353072 C 33040  01 
•0.371434242540361D  03 
•o.7oionoooonoooooD  03 
0.5716150328 8 120 30-14 


-0. 5751379273 C 75570- 
0.1387832077702^50 
n. 1013307718  <8  06  9. 3D 
”.2108°° 000° 000000 
-0. 157 r 4 42 3 3 10 00120- 
CT'-RUTED  GEOVALUE 


13  + 
03  + 
04  + 
03  + 
12  + 


>0.7783300030252500-14  I 
0.1410500057155310  03  I 
■0.  2700302337P.B25ED  03  I 
•0. 3330000000000000  03  ! 
■n.22n174700424B 040-13  I 


^ .  5  f>  ^  40  5  5  0  5 .8  0  5  8  0  B  0 

-p.  Tp-ipppppnpoin 
_  p  5 1 P  ^  0-  'i  n  P  ^  0  0  f)  p  0  0  0 

0  !  1  ''PPG  fif  f  C  GG  RG  f  70 
n. 00573^8:035500310 
RELATIVE  ERROR 
0.0000040223740010 
■8 . 4  5  7  0  C  C  9  9  7  6  5  7  8  7  ’  0  - 
0.4712447354 656920- 
9.4579009970578770- 
0.0 


1B  + 
P0  + 
00  + 
01  + 
00+ 


08 

15 

15 

15 


0 . 3  0  0  7  2  4  2  5  4  C 1  ?  0  f.  r>  D  1C  I 
0. 500OO0O0000O000D  00  i 
-0. 50^0000000000000  00  I 
-0.0666666656606600  00  I 
-0. 15270040591410 40  00  | 
RFLATIVE  RESIOUAL 
0. 4 0555 G 314306225 0-10 
1.°G17 73 5097136 4 50-10 
0.112 4. 3  91536 55 1190-15 
0. 1073640203077870-15 
•).  9209151010414300-16 


no.  OF  ITERATI-T 


59 


X 

II 


M 


9 


Example  2 


NONLINEAR  DIVISOR 
I  UP  MATRIX  A: 


311. + 

-330. 1 

307. + 

-003. 1 

OO  .  _ 

»**-»•  T 

1032.  1 

-250. + 

240.  1 

-112. +  -672 .  1 

64.+ 

-272.1 

-1053. + 

188 .  1 

127. + 

-604.  1 

?38.+- 

1000.  1 

183. +  -808  .  1 

-124. + 

420.  1 

-858.+ 

1023. 1 

81.  + 

176.  1 

670.  + 

-564.  1 

-3.+  1124.1 

00. + 

-776. 1 

531. + 

11.  1 

115. + 

-130.  1 

-377.+ 

843.  1 

-74.+  44.1 

-125. + 

-700. 1 

-SCO. + 

1057.  1 

153. + 

164.  1 

458.+ 

52.  1 

135. +  212.1 

THE  MATRIX  B: 

180. + 

11.1 

360. + 

47.  1 

72.+ 

104.1 

-180. + 

0.  1 

-144. +  96.1 

168.  + 

G.  1 

-546.+ 

GO.  1 

210. + 

46.  1 

588.+ 

24.  1 

-42.+  06.1 

00. + 

-61.1 

-450.+ 

-47.  1 

138. + 

-40.  1 

450. + 

-4.  1 

-6.+  -124.1 

42.+ 

20.  1 

336.+ 

-SI.  1 

-24.+ 

10.  1 

-258.+ 

-56.  1 

-60. +  4.1 

12. + 

77.* 

-354.+ 

-65.  1 

66.  ♦ 

-28.  1 

312. + 

-52.  1 

30. +  -28.1 

TRUE  EIGENVALUE 

1  n.iirpr»fin6fi6GGFfi7n  oi+ 

2  0.116RPGG6GGGF6G7D  01  + 

3  -o.ooooonoooooonnoo  oi> 

4  -o.onor.oonooonooneD  01+ 

5  -o.ORnononoonnoQoon  oi+ 
ALPHA 

1  -0.22022525938r724R  04+ 

2  -0 , G700835412374GGD  03  + 

3  “0.4228005411150180  03+ 

4  0.1800071442378100  03+ 

5  -0.014G81533504013D  01+ 

BETA 

1  “0. 1803072022-'.5820n  04  + 

2  -0.5527015071 4 50 f?D  03+ 

3  0.0575715471042510  02+ 

4  -0.2703234135710440  02+ 

5  0. 1310205083502520  00+ 
CLOUTED  EIGENVALUE 

1  0.1166666665837220  01  + 

2  o.iir.r.rr«rcr,74onnn  01+ 

3  -o.^ooooooooooooood  ni+ 

4  -".onooonoooooooooD  01+ 
3  -o.onooooooooooooon  01+ 

RELATIVE  ERROR 

1  .  7  n  3  3  3  3  7  5  f>  1 2  5  0  4  0  D  -  0  2 

2  0 . 7  03  5  3  8  7  2  0  4  7  0  HO  7  n-  0  8, 

3  0. 2 nl 227023 2133 100- 14 

4  0 . 5  3  0  2  4  5  0  4  5  7  2  3  0  8  2  D  - 1 5 

5  0 . 2  2  2  n  4  4  5  0  4  0  2  5  0  3 1  n  - 1 4 


-0.133333333333333D  01  I 
-0.133333333333333D  01  I 
0.1000000000000000  01  I 
-o.monnnoooooooooD  01  1 

-O.IOOOOOOOOOOOOOOD  01  I 

0.2531048717204740  04  I 
0.720*1040128040000  03  I 
0.1585880914052330  04  I 
0.  587',47003142127D  03  f 
0.7150311833064200  02  I 

0.5273030380225330  01  I 
-0.1344025824034030  02  I 
-0. 15OO1375135021 20  03  I 
-0. '■.252302707012030  02  I 
-0.7050440750S88030  01  1 

-0. 13353333500373FO  01  I 
-0.1333333310570310  01  I 
0. 000000000000080D  00  I 
-0.0000000900009070  00  I 
-0.00090099090097SO  00  I 

RELATIVE  RESIDUAL  NO.  OF  ITERATIONS 

0.53071104 6 2593 4 3 0-16  0 

0 .  5 550011 9 0 G P  *  4 2 0-1  *5  1 

0. IP 4 5715301074350-15  3 

2. 12745G207P05S55D-15  1 

o. 0010465227344600-10  1 


60 


9^s^l^«i^43QgBaSJ!3' 


SUBROUTINE  GELHES IND, N, A ,B , WAN7X ,X,EPSA, EPSBI 


THIS  SUBROUTINE  REDUCES  THE  COMPLEX  MATRIX  A  TO  UPPER 
H6SSENBERG  FOKH  ANO  REDUCES  THE  COMPLEX  MATRIX  B  TO 
TRIANGULAR  FORM 


C  INPUT  PARAMETERS: 

C 

C  ND  THE  RUW  DIMENSION  OF  THE  MATRICES  A,  B  «X 
C 

C  N  THE  OROER  OF  THE  PROBLEM 

C 

C  A  A  COMPLEX  MATRIX 

C 

C  BA  COMPLEX  MATRIX 

C 

C  WANT X  A  LUG1CAL  VARIABLE  WHICH  IS  SET  TO  .TRUE.  IF 
C  THE  EIGENVECTORS  ARE  WANTED.  OTHERWISE  IT  SHOULD 

C  BE  SET  TO  .FALSE. 

C 

C  OUTPUT  PARAMETERS: 

c 

C  A  ON  OUTPUT  A  IS  AN  UPPER  HESSENBERG  MATRIX,  THE 
C  ORIGINAL  MATRIX  HAS  BEEN  DESTROYED 

C 

C  fa  AN  UPPER  TRIANGULAR  MATRIX,  THE  ORIGINAL  MATRIX 
C  HAS  BEEN  DESTROYED 

r 

C  X  CONTAINS  THE  TRANSFORMATIONS  NEEDED  TO  COMPUTE 
C  THE  EIGENVECTORS  OF  THE  ORIGINAL  SYSTEM 

C 

C  EPSA  THE  NORM  OF  A*THE  PRECISION  OF  THE  MACHINE 
C 

C  EPSB  THE  NORM  UF  B*TH£  PRECISION  OF  THE  MACHINE 
C 

C****  THE  VALUE  OF  MACHEP  IS  MACHINE  DEPENDENT ******* 

(;*♦**♦!*  lS  SET  FUR  THE  IbM  360  MACHINE,  DOUBLE  PRECISION***** 
C 

C  PROBLEMS  WITH  THIS  SUBROUTINE  SHOULD  BE  DIRECTED  TO: 

C 

C  LINDA  KAUFMAN 

C  SERRA  HOUSE 

C  COMPUTER  SCIENCE  DEPARTMENT 

C  STANFORO  UNIVERSITY 

C 

COMPLEX  *16  Y, AIND, NOI , BIND, ND) ,XIND, NDI 
RcAL*o  ANi,BNI,C, RABStD, EPSA, EPSB, MACHEP, ANORM.BNOKM 
LOGICAL  WANTX 
NMl=N-l 
C 

C  COMPUTE  EPSA, EPSfa 


Av»*44.«aj«3i 


i  ' 

k  i 

a  i 

i  j 

? ! 
;  ! 
?  ! 
t  ! 

i  ; 

I  5 


c 

c 

c 


10 


11 


12 

15 


18 


19 


RABS ( A( I t I-l ) ) 


H AC  HEP *2. 220-16 
ANORM  *  0. 

BNORM  »  0. 

OU  5  1*1# N 
ANI  *  0. 

IF  (i.NE.l)  ANI 
BN1  *  0. 

00  3  J*l#N 

ANI  *  ANI  ♦  RABSU (I.JIi 
BN I  *  BNI  ♦  RABSIBI 1« J) I 
CONTINUE 

IF  (ANl.GT.ANORNI  ANORH  *  ANI 
IF  ( BN  I • GT •BNORM I  BNORM  *  BNI 
CONTINUE 

IF  (ANORM.EQ.O.I  ANORM  -  MACHEP  1 

IF  (BNORM.EQ.O.J  BNORM  *  MACHEP 
EPSA  -  MACH E P&ANORM 
EPSB  *  MACHEP*BNORM 

REOUCE  B  TO  TRIANGULAR  FORM  USING  ELEMENTARY  TRANSFORMATIONS 

00  30  1*1#NM1 

0*0.000 

IP1*UI 
00  10  X»IP1 »N 

C*KABS( B(K# II I 
IF  (C.LE.OI  GU  TO  10 

u*c 

II*K 

CONTINUE 

IF  (O.EQ.O.OOI  GO  TO  30 
IF  e O.LE.RABS ( B( 1(1)1)  GO  TO  15 
DO  l  J*l,N 
*  *A(l • J) 

A  1 1* J l*A( II, J) 
ftf  1 1 #  J  J*  Y 
00  12  J*I#N 

Y  -Btl » J  ) 

bu#jj*bui#j) 

B(  lit J)*Y 
00  20  J*IP1 » N 

Y*B(J#I)/B(i#I) 

IF  (RABS  t Y) .EU.O.OO)  GO  TO  20 
00  18  K*1 »N 

A(  J  #K)*A( J  »K|— Y*A( I »  K1 
00  19  K*  IPl »N 

B( J.X)»B(J,K)-V*B(I,K> 

8(J»I »=(0.00#0.00) 

CONTINUE 

CONTINUE 


20 
30 
C 

C  INITIALIZE  X 
C 


IF  (.NOT.WANTX)  GO  TO  AO 


62 


T-,izy- 


37 

38 
C 
C 

c 

40 


45 


46 

50 


52 


54 

C 

58 


60 


64 


68 

70 


72 


74 


00  38  1=1, N 

00  37  J*i,N 

X(I,Ji  =  (0.DO,r,.O0i 
XCI, 11=11.00,0.000) 

REDUCE  A  TO  UPPER  HESSENBERG  FORM 

NM2*N-2 

IF  CNM2.LT. W  GO  TO  100 
00  90  J=1 »NM2 
JM2*NM1- J 
JPi=J«-i 

00  80  l 1*1 , JM2 
I=N+1-II 
IM1=I-1 

IF  (RABSC  AC  it J J I.LE.RABSC AC IMI ,Ji ) I  GO  TO  50 
DO  45  K*J,N 

Y=AC I ,K I 

AC  1 ,K)=A( IMI »K) 

A(IM1»K)=Y 
00  46  K=IM1,N 
Y=B(  I.»K) 

8(1,  K  )=B  C  I  Ml  ,K.) 

8C lMi, K1=Y 

IF  (RABS(Ali,J)  I.EQ.O.OOi  GO  TO  58 
Y=A(I»JI/A( IH1 , J ) 

00  52  K*JPi,N 

AC  i  ,K)=ACI  ,K)-Y*AUM1,K) 

All, J)=CO. 000, 0.001 
00  54  K=1 Ml , N 

BCI  ,K)=B(I,K)-Y*BUMl,Kf 
TRANSFORMATION  FROM  THE  RIGHT 

IF  CRABSC B( 1 , 1M1) ) «LE .RABSC 6(1,1)))  GO  TO  70 
00  60  K=l* I 

Y*6C  K, 1 1 

B( K» l J  =B( K  »  IMi ) 

B ( K, I  Ml ) =Y 
DO  64  K=1,N 

Y=ACK,I) 

AC  K, I ) *A( K  » IMI  ) 

ACK,IM1)*Y 

IF  ( .NOT .WANT X)  GO  TO  70 
00  68  K=1,N 

Y=X(K, n 
X  CK ,  1  )=XCK,IM1) 

X ( K , I  Mi) *Y 

IF  CRABS CBC I, CHlii.EQ.O.DO)  GO  TO  80 
Y=8CI ,IMl)/B(I,l) 

00  72  K= 1 , IMI 

BCK, IMi)=B(K, IMI )-Y*B( K,  I ) 

B  (  I  ,IMU  =  C0. 00,  0.001 
OU  74  K=1,N 

A(K,  IM1)-ACK,IM1 )-Y*A C K, I ) 

IF  i .NO  T. WANTX )  GO  TO  80 
00  76  K= 1 , N 


6;» 


76  X(K,IMl)*A<fc»lMD-V*X<K*I  ) 

80  CONTINUE 

90  CONTINUE 

100  RETURN  1 

ENO 

MJBRUUlINE  GLR(ND*NN* A,B***X» ITER,EPSA*EPS8*MANTX»EiGA,EI6B) 

C 

C  THIS  SUBROUTINE  SOLVES  THE  GENERALIZED  EIGENVALUE  PROBLEM 
C  AX*  LAMBDA  B  X 

C  WHERE  A  IS  A  COMPLEX  UPPER  HESSENBERG  MATRIX  OF  ORDER  NN  AND  B  IS 
C  A  COMPLEX  UPPER  TRIANGULAR  MATRIX  OF  ORDER  NN 
C 
C 

C  INPUT  PARAMETERS 
C 

C  | 

C  ND  ROW  DIMENSION  OF  THE  MATRICES  A ,B*X, I  TER ,EIGA»EIGB 

C  NN  ORDER  OF  THE  PROBLEM 
C 

C  A  AN  NN  X  NN  UPPER  HESSENBERG  COMPLEX  MATRIX 
C 

C  B  AN  NN  X  NN  UPPER  TRIANGULAR  COMPLEX  MATRIX 
C 

C  *  ERROR  RETURN?  IF  AFTER  iO  ITERATIONS*  THE  NORM  OF  THE 
C  SUBDIAGONAL  UF  A  HAS  NOT  SHOWN  A  SUFFICIENT  DECREASE 
C 

C  X  CONTAINS  TRANSFORMATIONS  TO  OBTAIN  EIGENVECTORS  OF 
C  ORIGINAL  SYSTEM 

C  IF  GELHES  HAS  NOT  BEC.N  USED*  X  SHOULD  BE  THE  IDENTITY  MATRIX 
C 

C  WANTX  LOGICAL  VARIABLE  WHICH  SHOULD  BE  SET  TO  .TRUE.  IF  EIGENVECTORS 
C  ARE  WANTED.  OTHERWISE  IT  SHOULD  BE  SET  TO  FALSE 
C 

C  EPSA  THE  NORM  UF  A  TIMES  THE  MACHINE  PRECISION.  NEED  NOT  BE 
C  SET  IF  GELHES  HAS  BEEN  USED 

C  EPSB  THE  NORM  OF  B  TIMES  THE  MACHINE  PRECISION.  NEED  NOT 
C  BE  SET  IF  GELHES  HAS  BEEN  USEO 

C 
C 

C  OUTPUT  PARAMETERS 

C 

C 

C  W  THE  ITH  COLUMN  CONTAINS  THE  ITH  EIGENVECTOR  IF  EIGENVECTORS  ARE 
C  REQUESTED 

C 

C  ITER  AN  INTEGER  ARRAY  OF  LENGTH  NN  WHOSE  ITH  ENTRY  CONTAINS  THE  NUMBER 
C  OF  ITERATIONS  NEEDED  TO  FIND  THE  ITH  EIGENVALUE 

C 

C  EIGA  AN  NN  ARRAY  CONTAINING  THE  DIAGONAL  OF  A 
C 

C  E1G8  AN  NN  ARRAY  CONTAINING  THE  DIAGONAL  OF  B 
C 

C  THE  ITH  EIGENVALUE  CAN  BE  FOUNO  BY  DIVIDING  EIGA1I)  BY  EIGBIII 
C  WATCH  OUT  FOR  EIGBII)  BEING  ZERO 
C 

C***«**«*  THE  QUANTITY  MACHEP  IS  MACHINE  DEPENDENT  ********* 
t*********  it  is  SET  FOR  THE  IBM  360,  DOUBLE  PRECISION******** 


PROBLEMS  WITH  THIS  SUBROUTINE  SHOULD  BE  DIRECTED  TO 
LINDA  KAUFMAN 

SERkA  HOUSE,  COMPUTER  SCIENCE  DEPARTMENT 
STANFURO  UNIVERSITY 


L0MPLEX*i6  A(ND,NOi ,B (NO, NO) , EIGA(ND) ,E IGBINDI 
COMPLEX* 16  S,T,W,Y,Z,OCMPLX,CDSGRT,TS 
INTEGER  IT  ER(NU) 

C0MPLEX*L6  ALFM,BETM«O«SL,0ENt  ANN«ANM1N, ANM1M1 
REAL*8  EPSA ,EPSB, SS ,R ,OLD ,NEW 
CGMPLEX*16  X ( NO, NO  I 

RE AL* 8  MACHE P / 2. 220-1 6/ , 00 , 01 ,02 • EO , E l ,R ASS , OABS 

LOGICAL  WAN TX 

N=NN 


IF  IN.LE.l)  GO  TO  10Q 
I  T  S=0 
NM1=N-1 
00  12  L B=2, N 
L*N*2-LB 

IF  I RABS (AIL, L- II ) »LE .MACHEP*! RABSi A(L— 1 »L— 1 ) I 
K  +RABS(A(L,L) J))  GO  TO 

CONTINUE 
L=1 

IF(L.EW.N)  GO  TO  100 
IF  (ITS.LT.30)  GO  TO  20 
IF  (ITS.GT.jOJ  GO  TO  16 
OLO=O.DO 
DO  15  1=1, NM1 

OLD=OLD-*-RABS  I  A(  1*1,11) 

GO  TO  20 
NEW=Q.00 
00  19  1=1, NMi 

NEW=NEW+RABS( Al 1*1,11 ) 

IF  CNEW.GT.0.5*0L0I  RETURN  1 
OLD=NEW 


C 

C  CHECK 
C 


CONSECUTIVE  SMALL  SUBDIAGONAL  ELEMENTS 


IFCN.EU.L  +  U  GO  TO  25 
02= RABS (A{N-i,N-i! 1 
£ 1=RABS(A(N»N— 1 ) ) 

01= RABS (AIN, Nil 
NL=N-(L*1I 
DO  24  MB=1,NL 
M=N-M8 
E0=E  1 

E  l=RABS(  A(M,M-U  ) 

00=01 

01=02 

02=RA8Sl AIM-1, M-l) I 
IFCEO*£l.LE.MACHEP*0l*I00*0l*02) ) 
CONTINUE 


6' 


2  b  M*L 

26  CONTINUE 
C 

IFl ITS.EQ.10.0R.ITS.EU.20)  GO  TO  38 
C 

C  COMPUTE  SHIFT  AS  EIGENVALUE  OF  LOWER  2  BY  2 
C 

ANN*A(N»N) 

ANMIN=A|NHI,N) 

ANMIM1*A(NM1,NMI) 

S= ANN*8(NM1 »NM1 )-(A(N»NMl))*8(NMl»N) 
ta*A(N,NMl)*BlN,N)*(ANMIN*B(NMlfNMl)- 
1  B( NMl»  N)*ANM1M1 } 

Y=(ANMiMl*B(N,N)-S)/2. 

ZsCOSQKT I Y* Y*W) 

IF  (RABSCZ) .EQ.O.DO)  GO  TO  36 
DQ*Y/Z 

IFIOQ.LT. 0.00)  Z*-Z 
36  DEN*(Y*Z)*B I NM1 »NM1 N) 

W=A(M,M)*DEN-B(M,K)*( IY*Z)*S-W) 
ZxA(M+l,M)*OEN 
GO  TO  40 
C 

C  AO-HOC  ShIFT 
C 

38  W*A4  N» N— l ) 

Y=*A4  N-l  *N-2  ) 

W*A4N» M|— DCNPLX4RABS4  W) »RA8St Y) ) *BC M,M) 
Z=A(M*1,M) 

40  CONTINUE 
ITS-  ITS+l 
C 

C  FIND  L  ANO  M  ANO  SET  A3 LAM  ANO  B*LBM 
C 

NP l*N+l 

LURI=L 

NNQRN*N 

IF  ( .NOT .WANTX)  GO  TO  42 
LORI3! 

NNURN=NN 

42  00  9C  I  =  Mf  NM1 

J*I  +  l 
C 

C  FIND  ROW  TRANSFORMATIONS  TO  RESTORE  A  TO 
C  UPPER  HESSENBERG  FORM.  APPLY  TRANSFORMATIONS 
C  TO  A  ANO  B 
C 

IF  II.EU.M)  GO  TO  bO 
M-A4 1 » I-iJ 
Z*A4J# i-1) 

bO  IF  (RABS(W).GE .RABS ( Z) )  GO  TO  6G 

C 

C  MUST  PIVOT 
C 


Du  55  K*I t NNORN 


Y=A( I,K) 

A(1 |K)*A( J  ,K1 
A ( J , K)= Y 
V=B( I ,K) 

BU  ,KI*BIJ»K) 

55  BiJ,K)=Y 

Y*W/Z 

if  tl.GT.M)  Ai I,I-l)*AtJ,l-ll 
GO  TO  62 
60  V=Z/K 

62  if iRABStY) .EU.O.DQ)  GO  TO  65 

00  64  K=I,NN0RN 

Ai J  ,Ki=AiJ,K)-Y*AiI,K) 

64  BiJ  ,K)=6tJ,K)— Y*6ii»K) 

if  il.GT.Mj  Ai  J,I-l)=iO.DO,O.DO) 

C 

C  PERFORM  TRANSFORMATIONS  FROM  RIGHT  TO  RESTORE  B  TO 
C  TR IANGLULAR  FORM 
G  APPLY  TRANSFORMATIONS  TO  A 
C 

65  If  tRABSiBi Jtiii.EO.O.OO)  GO  TO  11 

If  iRABSIBiJ*JJ) .GE.KABS iBi J, till  GO  TO  80 
C 

C  MUST  PIVOT  COLUMNS 
C 

00  70  K-LORUJ 
Y=A IK,  J  J 
AIK  »  J) =At  K» I ) 

A (K  » I  )  =  Y 
Y=B ( K, J  J 
BIK, J)=B{K,i) 

70  BiK  t  I  J=Y 

If  U.EQ.NMl)  GO  TO  75 
Y=At J*l, Ji 
A( J+i, J»  =  .4t J*i, I) 

Ai J+i,I >=Y 

75  IFt .NUT.WANTX)  GO  TO  BO 

DO  7ti  K= l »  NN 
Y=XIK, Ji 
XiK, J)=XIK,  I ) 

78  X( K , I i=Y 

80  If  iRABSiBU.iii.EU.O.OO)  GO  TO  90 

Z=8( J , l  )/Bi J  *  J ) 

00  82  K=LOKi,J 

AiK»I)=AiK»Ii— Z*AtK»J) 

82  BiK,I)=BlK,I)-Z*8(K,J) 

BiJ, I  i  =  10.00,0.00) 

If  ti.LT.NMl)  Atd2,I)*AiU2,l)-Z*Ail»2,J) 
I  f  (  •  NO  T  •  HAN  IX  )  GO  TO  90 
00  «5  K=i,NN 

85  X(K,1 j=XiK, I)-Z*X<K,J) 

90  CONTINUE 

GO  TO  11 
C 

100  CONTINUE 


i 

t 


6? 


EIGA(N)=A(N,N) 

EIGB(N)=B(N,N) 

IF  (N.EO.i)  GO  TO  HU 

ITERINl^ITS 

N»NM1 

IF  (N.GT.L)  GU  TO  10 
ITER( 1)*Q 
GO  TO  100 
C 

C  FIND  EIGENVECTORS  USING  B  FOR  INTERMEDIATE  STORAGE 

C 

110  IF(.NOT.WANTX)  RETURN 
M»NN 

115  CONTINUE 

Al  rM*A(M,M) 

BEfM*B(M,M) 

B(M,M)=(l.OOtO.OO> 

L  *  M-l 

IF  (L.EU.O)  GO  TO  140 

120  CONTINUE 

LI  *  L*1 
SL  *  0. 

00  130  J-LlfM 

SL  =  SL  ♦  (BET M*A ( L *  J ) -ALF M*  B(  L  »  J) )  *B  C  J  »  M) 
130  CONTINUE 

0  =  BETM*A(L»L|-ALFM*B1L, L) 

IF  (RABSIDI.EQ.O.)  0  «  (EPSA+EPSB) /2. 

B(L,Mi  *  -SL/O 
L  =  L-l 

140  IF  (L.GT.O)  GO  TO  120 

M=M-1 

IF  (M.GT.O)  GO  TO  115 
C 

C  TRANSFORM  TO  ORIGINAL  COORDINATE  SYSTEM 
C 

M  =  NN 

200  CONTINUE 

00  220  1 =1»NN 
S  *  0. 

00  210  0=1  #M 

S  =  S  ♦  X(I,J)*B(J,M) 

210  CONTINUE 

X ( I ,M  )  *  S 
220  CONTINUE 

M  =  M-l 

IF  (M.GT.O)  GO  TO  200 
C 

C  NORMALIZE  SO  THAT  LARGEST  COMPONENT  *  1. 

C 

M  *  NN 

230  CONTINUE 

SS  =  0. 

00  235  1=1, NN 

R  =  KABS( X( I ,M)  ) 

IF  (R.LT.SS)  GO  TO  235 


68 


SS  *  R 
D  *  XU, Mi 
CONTINUE 

IF  (SS.EQ.O.UOI  GO  TO  245 
DO  240  I=l,NN 

XiI,M)  *  X( I »  Mi/D 
CONTINUE 
M  =  M-l 

IF  (M.GT.Oi  GO  TO  230 

RETURN 

END 


REAL  FUNCTION  RA3S*8(Z) 
C0MPLEX*16  I.i  ll 
REAL*8  Tt2),0ABS 
EQUIVALENCE  (ZZ.TUi) 
zz*z 

RAttSs0A8S(T I  li i*DA8S( Tl  21 ) 

RETURN 

ENO 


i  A  '■  :.»<r/^AlW  A;  r^rr/^kV-A^ybLyV^ti  ? n 


J 


4 


■t— 


4 


4 


4 


4 


Acknowledgements 

I  am  deeply  indebted  to  my  advisor  Professor  Cleve  Moler 
for  his  guidance,  encouragement,  and  enthusiasm.  I  thank  him  for 
patiently  explaining  to  me  the  intricasies  of  the  QZ,QP,  pr..\  ip 

i 

algorithms  and  for  challenging  me  to  find  the  LZ  algorithm.  I 
am  grateful  to  Professor  George ' Forsythe ,  Professor  John  Herriot,and 
Professor  Gene  Golub  for  reading  the  manuscript  and  to  Professor 
Michael  Osborne  for  finding  an  error  in  one  of  the  earlier  versions 
of  one  of  the  convergence  proofs.  I  would  also  like  to  thank 
Mary  Bodley  and  Fran  Brumbaugh  for  typing  the  manuscript. 


70 


/ 


REFERENCES 

1.  H.J.  Buurena,  "A  Geometric  Proof  of  the  Conver' -'•nee 

of  the  OP  Method"/ Report  TW-62,  Metl-eniCtl ’’  .»i 
Insti  tuut, Groningen, North  East  Netherlan  s,1968. 

2.  C.R.  Crawford, "The  Numerical  Solution  of  toe 

Generalized  Eigenvalue  Problem"/  University 
of  Michigan  Ph.D.  Thesis, 1971. 

3.  G.  Fix  and  R.  He! herder, "An  Algorithm  for  the 

Ill-conditioned  Generalized  Eigenvalue 
Problem"/ to  be  published  In  Num.  Math. 

4.  G. Forsythe  and  C.  Mol er. Computer  Solution  of  Linear 

A1 eebralc  Systems, Prentlce-Hal 1 . Inc. /Englewood  Cliffs, 
New  Jersey/ 1967. 

5.  J.G.  Franc!S/"The  0.R  Transfornatlon-A  Unttary 

Analogue  to  the  LR  Transformation"/ 

Computer  Journal  4//265-271, 332-345/1961-1962. 

6.  F. R.Gantnacher . App.l.i cation  of  the  Theory  of  Matrices. 

Intersciences  Publ i shers, I nc. ,  New  York/ .959. 

7.  G.H.  Golub/R. Underwood/ J.H.  Wil kl nson, "The  Lanczos 

Algorithm  for  the  Synnetrlc  Ax^Bx  Problem"/ 

Stanford  Computer  Science  Report  270/March  1972. 

8.  P.  Lancaster ,  JjaobJjl  Mgtr  j.£.g,S.  ami  V.,1  bCdt  LilX  SYS.tfimS/ 

Pergamon  Press, New  York, 1966. 

9.  R.S.  Martin  and  J.M.  Wi  1 ki nson, "Reduct  Ion  of  the 

Symmetric  Ei genprobl em  Ax  = A  Bx  and  Related 
Problems  to  Standard  Form". Numer .  Math.  11, 

99-110/ 19F3. 

10.  C.B.Moler  and  G.W.  Stewart, "An  Algorithm  for  the 

Generalized  Matrix  Eigenvalue  Problem  Ax»ABx", 

Stanford  Computer  Science  Report  232, August  1971. 

11.  B.N.  Parlett/"GlobaI  Convergence  of  the  Basic  OR 

Algorithm  for  Hessenberg  Matr I ces".Math.  Como.  22, 
803-317/1968. 


— and  W,  Kahan,"0n  the  Convegrence  of  the  Practical 
OR  Algorithm"  .Proceedings  of  the  I F  IP  Congress  196.8. 
A25-A30 . , 1968 . 


— and  W.G. Poole, "A  Geometric  Theory  for  the  QR, 

LU,  and  Power  I  terat  l ons", Berkel ey  Computer  Science 
Technical  Report, 1971. 


14.  B.N.  Pari ett,  "Canon leal  Deconposl  t  Ion  of  Hessenberg 

Matrices”. Math  Comp.  21,223-227,1967. 

15.  G. Peters  and  J.H.  Wt 1 kl nson, "Ax  =^Bx  and  the  Generalized 

Eigenorohl em",Slam.  J.  of  Muner.  Anal.  7,479-492,1970. 

16.  G.  Peters  and  J.H.  Wilkinson, "Eigenvectors  of  real  and 

complex  matrices  by  LR  and  QR  trlangularlzatl^ns", 
Numer.  Mathu  16,191-204,1970. 

17.  H.  Rutishauser,  "Solution  of  the  Eigenvalue  Problem  with 

the  LR  Transformation",  National  Bureau  of  Standards 
Applied  Math  Series  49,  January  1958. 

18.  G. W. Stewart, "On  the  Sensitivity  of  the  Eigenvalue  Problem 

Ax=  ABx",  Center  of  Numerical  Analysis,  University  of 
Texas  Report  13,  March  1971. 


72 


